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CHAPTER  1 ; REVIEW  OF  CLASSICAL  FOURIER 


From  the  mathematical  point  of  view,  the  key  concept  in  Fourier 
methods  of  analysis  is  that  a more  or  less  arbitrary  function  can  be 
expressed  as  a (generalised)  sum  of  the  trigonometric  functions  (sines 
and  cosines)  or  their  equivalent,  the  complex  exponentials.  The  sum 
may  take  the  form  of  an  infinite  series  (classical  Fourier  series) , an 
integral  (classical  Fourier  transform) , or  a finite  series  (discrete 
Fourier  transform) . Classical  Fourier  methods  concentrate  on  the  case 
where  the  arbitrary  function  is  given  by  some  analytical  expression  whilst 
present  interest  is  more  with  a function  expressed  as  a discrete  set  of 
numerical  values;  the  latter  format,  of  course,  is  especially  appropriate 
for  computer  applications.  However,  before  discussing  what  can  perhaps 
be  best  termed  "numerical  Fourier  methods"  which  are  based  on  the  discrete 
Fourier  transform  it  is  useful  to  review,  albeit  briefly,  their  classical 
counterparts.  Of  the  very  extensive  literature  on  the  classical  methods 
special  attention  is  drawn  to  the  accounts  by  Hsu(l),  Champeney (2) , 
Bracewell(3)  and  Papoulis(4);  these  are  all  oriented  towards  applications 
of  current  interest.  Other  standard  texts  are  listed  in  references  (5),  (6) 


1.2  Physical  Motivation  for  Fourier  Series 


The  fact  that  an  arbitrary  function  can  be  expressed  as  a sum  of  the 
trigonometric  functions  can  hardly  be  deemed  intuitively  evident  and  it  may 
be  useful  to  describe  briefly  how  such  a result  was  discovered.  Essentially 
it  came  out  of  studies  of  the  vibration  of  a tense  string  carried  out  by 
Daniel  Bernoulli  around  1750.  (Truesdell(9)  has  described  at  length  the 
"celebrated  and  deplorable  quarrel"  which  arose  between  Bernoulli  on  the  one 
hand  and  Euler,  Lagrange  and  D'Alembert  on  the  other,  over  the  validity  of 
this  work) . 


The  partial  differential  equation  governing  the  small  vibrations  of 
a uniform  string  of. line  density,  d,  and  under  a tension,  T,  is 
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where  u denote  the  lateral  displacement  at  a distance  x along  the  string 
2 

at  time  t,  and  c = T/d.  Confining  attention  to  the  case  when  the  string 
is  initially  at  test  in  some  deformed  configuration  there  will  be  initial 
conditions  of  the  form 

u(x,o)  - g(x)  (1.2.2) 

ut(x,o)  = 0 (1.2.3) 

As  well,  there  will  be  boundary  conditions  at  the  ends  x=0  and  x=L  of  the 
string.  First,  consider  the  case  of  ’fixed  ends'  Fig. 1.1  where 


u(o,t)  * u(L,t)  = 0 (1.2.4) 

Use  of  the  method  of  "normal  modes"  (which  is  here  equivalent  to  the  method 
of  "separation  of  variables" ) shows  that  an  expression  of  the  form 
00 

u(x,t)  » Z b sin  (nJIx/L)  cos  (nJIct/I*)  (1.2.5) 

i n 
n=l 

satisfies  the  conditions  (1.2.1),  (1.2.3)  and  (1.2.4).  (In  fact,  each  term 
in  the  series  satisfies  these  conditions).  The  condition  (1.2.2)  now  becomes 
00 

g(x)  = E b sin  (nllx/L)  ; (1.2.6) 

n=l  n 

but,  since  g is  more  or  less  arbitrary,  eqn  (1. 2£)  implies  that  an  arbitrary 
function  can  be  expressed  as  an  infinite  series  of  sines.  Assuming  the 
validity  of  the  expansion,  the  coefficients  bn  can  be  readily  determined  by 
using  the  "orthogonality  relations",  namely,  that  for  integral  m,n 

L 

/ sin  (mllx/L)  sin  (nllx/L)  dx  = o,  m^n  (1.2.7) 

o 

= L/2,  m=n  (1.2.8) 

(This  is  an  example  of  the  general  orthogonality  of  normal  modes  of  vibration) 
Hence, 

b - (2/L)  /L  g(x)  sin  (nirx/L)  dx  (1.2.9) 


and  the  problem  is  now  solved. 

) 

One  other  result  might  be  noted.  On  multiplying  each  side  of  eqn 
(1.2.6)  by  itself  and  integrating,  it  can  be  established  that 


f L{g(x)  ) 2 dx 
o 


(t/2)  Z 

n=l 


n 


(1.2.10) 
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This  result  is  generally  referred  to  as  Parseval's  theorem  for  a sine 
series. 

Now  consider  another  set  of  boundary  conditions,  in  particular 
the  case  of  "free  ends"  where,  instead  of  eqns  (1.2.4),  we  have 


ux  (o, t)  = ux(L,t)  = 0 


(1.2.11) 


' (This  slightly  artificial  case  can  be  motivated  by  considering 
the  ends  of  the  string  as  attached  to  springs  (Fig.  1.2)  when  the  boundary 
conditions  are  of  the  form 


T u - Ku  = 0 
x 


(1.2.12) 


where  K is  the  stiffness  of  the  springs;  for  vanishingly  small  K,  eqns 
(1.2.11)  result.)  Separation  of  variables  now  leads  to  an  expression 
of  the  form 


u(x,t)  = a /2  + E a cos(nIIx/L)  cos(nIIct/L) 
o , n 

n=l 


(1.2.13) 


This  satisfies  (1.2.1),  (1.2.3)  and  (1.2.11).  The  condition  (1.2.2)  becomes 


g(x)  = a /2  + E a cos(nIIx/L) 
n=l 


(1.2.14) 


(The  point  of  associating  the  factor  of  1/2  with  the  constant  term  - which 
physically  corresponds  to  a rigid  body  motion  - will  appear  shortly.)  Thus, 
it  appears  that  an  arbitrary  function  can  also  be  expressed  as  an  infinite 
series  of  cosines.  The  "orthogonality  relations"  are  now  that,  for 
integral  m,n  (including  zero) 


/L  cos  (mllx/L)  cos(nIIx/L)  dx  )=  0 , m / n 
o 

= L/2  , m = n / 0 


= L , m = n = 0 


Hence, 


a = (2/L)  f g(x)  cos  (nllx/L)  dx 
n o 


(1.2.15) 

(1.2.16) 


(1.2.17) 


(1.2.18) 


(Use  of  the  factor  of  1/2  in  (1.2.14)  means  that  eqn  (1.2.18)  can  be  used 
for  all  integral  values  of  n including  zero.) 

Parseval's  theorem  for  a cosine  series  is  as  follows  : 


/L  {g(x)}2  dx  = (L/2)  (aV2  + E ah 


(1.2.19) 


As  already  remarked  the  above  arguments,  for  the  sine  series  at 
any  rate,  are  essentially  due  to  Bernoulli.  Fourier's  contribution 
(10)  made  in  1822  consisted  in  the  elucidation  of  many  points  left 
obscure  by  Bernoulli,  and  included  the  actual  derivation  of  numerous 
series;  also,  he  was  the  sole  discoverer  of  the  Fourier  integral  to  be 
discussed  later. 

1.3  The  Different  Fourier  Series 

1.3.1  Fourier  Sine  Series 

Recalling  the  argument  for  a string  with  fixed  ends,  the 
assertion  is  that  a fairly  general  function  g(x)  defined  over  the  interval 
(0,L)  can  be  represented  in  the  form 

00 

g(x)  = Z b sin  (nllx/L)  (1.3.1) 

n=l  n 

where  b = (2/L)  / *g(x)  sin  (nJIx/L)  dx  (1.3.2) 

n o 

As  an  example,  consider  the  function 

g(x)  = exp  (x/L)  - 1 (1.3.3) 

which  is  shown  in  Fig. 1.3 (a).  Here  the  Fourier  coefficients  are 

b = 2 {(-l)n+1  { (e-1)  n2n2  -i}  -1  } / {nil  (l+n2II2)}  (1.3.4) 

n 

The  representation  of  the  function  that  is  recovered  by  summing  (1.3.1) 
up  to  n=4  is  shown  in  Fig  1.3(b).  It  might  be  noted  that  eqn  (1.3.3) 
violates  one  of  the  conditions  of  the  string  problem  inasmuch  as  g(L)  is 
non  - zero.  Whilst  the  series  representation  with  only  four  terms 
reproduces  the  function  at  least  fairly  over  most  of  the  range  it  is  badly 
in  error  near  x = L;  there  it  actually  returns  the  value  zero  instead  of 
(e-1). 


Although  in  the  string  problem  it  is  of  no  concern  how 
the  Fourier  series  behaves  outside  the  interval  (0,I>),  nevertheless  it  is 
informative  to  look  at  this.  Firstly,  since  the  sine  is  am  odd  function 


the  series  (1.3.1)  always  returns  the  result 


g(-x)  = -g(x)  (1.3.5) 

Secondly,  since  for  any  integer  n, 

sin  {nil  (x  + 2L)A>  = sin  (nllx/L)  (1.3.6) 

the  series  returns  the  result 

g(x  + 2L)  = g(x)  (1.3.7) 

i.e.,  the  sine  series  represents  a periodic  function  of  period  2L.  Hence, 
for  the  particular  function  (1.3.3)  the  series  will  really  be  representing 
the  extended  function  of  Fig.  1.3(b).  This  latter  has  jump  discontinuities 
at  x = L,  -L  , 3L,  -3L  etc.  and,  in  returning  the  value  zero  there,  the 
series  is  actually  taking  the  mean  of  the  values  on  either  side  of  the 
discontinuity.  Also,  in  the  neighbourhood  of  the  discontinuity  the  series 
with  a finite  number  of  terms  oscillates  quite  widely  about  the  true 
function;  this  effect,  known  as  the  Gibbs  phenomenon,  is  discernible  in 
Fig.  1.3(b). 

1.3.2  Fourier  Cosine  Series 

Similarly,  the  discussion  of  a string  with  free  ends 
suggests  that  a fairly  general  function  g(x)  defined  over  (0,L)  can  be 
represented  in  the  alternative  form 

oo 

q (x)  = a /2  + E a cos  (nllx/L)  (1.3.8) 

n=l  n 

where 

a = (2/L)  I L g(x)  cos  (nllx/L)  dx  (1.3.9) 

n o 

This  cosine  series  can  be  just  as  well  used  as  the  previous 
sine  series.  On  applying  it  to  the  function  (1.3.3)  it  is  found  that 

a = 2 (e-2) , a = 2 { (-l)ne-l}  / (1  + n2It2)  (1.3.10) 

o n 

The  representation  that  is  recovered  by  summing  the  cosine  series  up  to 
n = 4 is  shown  in  Fig.  1.3(c). 


The  representation  given  outside  the  basic  interval 
(0,L)  by  the  cosine  series  is  determined,  firstly,  by  the  fact  that  the 
cosine  is  an  even  function  so  the  series  always  returns  the  result 

g(-x)  = g(x)  (1.3.11) 

Secondly,  the  series  has  a period  of  2L  so  that  (1.3.7)  again  applies. 

Thus,  for  the  function  (1.3.3)  the  series  is  really  representing  the 

extended  function  of  Fig.  1.3(c).  In  this  case  it  happens  that  the 

extended  function  is  continuous  so  that  summing  to  only  n = 4 gives  quite 

a good  representation,  the  coefficients  a^  decreasing  more  rapidly  than 

the  b . However,  it  must  be  observed  that  the  zero  slope  end  condition 
k 

imposed  in  the  string  problem  does  not  apply  to  (1.3.3)  at  either  x = o 
or  x = L.  The  condition  is,  in  a sense,  approximately  satisfied  at  x = o 
because  of  the  even  character  of  the  extended  function. 

1.3.3  Fourier  Trigonometric  Series 

If  one  is  only  concerned  with  the  representation  of  a 
function  over  the  single  interval  (0,  L)  — and  this  is  the  situation  for 
many  boundary  value  problems  - then  either  the  sine  or  cosine  series  may 
be  used  as  convenient.  However,  if  one  is  concerned  with  some  sort  of 
periodic  phenomenon  over  the  whole  range  of  x neither  representation  may 
suffice.  For  example,  consider  again  the  function  given  by  eqn  (1.3.3)  and 
suppose  that  it  is  required  to  reproduce  this  function  over  each  interval 
of  length  L for  all  x i.e.,  the  extended  function  of  Fig.  1.3(d)  is  to  be 
represented.  As  already  seen,  neither  the  sine  nor  cosine  series  can  do 
this;  however  it  is  possible  to  develop  a combined  sine  and  cosine  series 
which  will.  The  argument  which,  incidentally,  is  due  to  Fourier  is  as 
follows. 

Any  function  g (x)  can  be  written  in  the  form 
g(x)  * E(x)  + 0(x)  (1.3.12) 

where 


E(x)  = (1/2)  {g(x)  + g(-x)}  (1.3.13) 

0(x)  **  (1/2)  (g(x)  - g(-x) } (1.3.14) 


! 


Clearly,  E(x)  is  an  even  function  and  0(x)  an  odd  function.  In  the  present 
case  g(x)  is  required  to  be  periodic,  with  period  L i.e.. 


It  follows  that  E (x)  and  O(x)  are  each  also  periodic  with  period  L 


For  the  moment  consider  just  E (x) . It  is  an  even  function 
hence,  from  the  considerations  of  Section  1.3.2  above  can 


E(x)  - a /2  + la  cos(2IInx/L) 
° n-1  n 


E(x)  cos  (2nJIx/L)  dx 


(Note  that  it  is  necessary  to  replace  the  "L"  of  eqns  (1.3.8)  and  (1.3.9) 
by  "L/2”  here  since  those  equations  related  to  a function  of  period  2L,  not 
L as  is  required.)  Substituting  from  eqn  (1.3.13)  into  eqn  (1.3.17),  and 
using  the  periodicity  relation,  leads  to  the  result 


(2/L) / g (x)  cos  (2nHx/L)  dx 


In  a similar  fashion  it  can  be  shown  that 


0(x)  = E b sin  (2nIIx/L) 

i n 

n=l 


b = (2/L)  I g (x)  sin  (2nIIx/L)  dx 


Summing  up  then,  a more  or  less  arbitrary  function  of 
period  L can  be  represented  in  the  form 


g(x)  = a /2  + E { a cos(2nIIx/L)  + b sin  (2nIlx/L)}  (1.3.21) 
o , n n 

n=l 


where  a and  b cure  given  by  (1.3.1%)  and  (1.3.20)  respectively.  Assuming 
n n 

the  validity  of  the  expansion  (1.3.21),  the  coefficients  a^  and  bR  can  also 
be  determined  directly  using  the  orthogonality  relations  for  the  combined 
cosine  and  sine  function  set,  namely. 
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where  u and  v denote  cos(2mrx/L)  and  sin  (2mrx/t)  respectively, 
n n 

Parseval’ s theorem  takes  the  form 


/ L{g(x) }2  dx  = (h/2)  {a2/2  + l (a2  + b2  ) } 
° ° n=l  n n 


(1.3.23) 


The  representation  (1.3.21)  is  the  general  form  of  a 
Fourier  series;  the  cosine  and  sine  series  can  be  regarded  as  particular 
cases  appropriate  to  even  or  odd  functions  respectively. 

The  series  (1.3.21)  could  have  been  derived  from  physical 
considerations  by  solving  eqn  (1.2.1)  with  initial  conditions  (1.2.2)  and 
(1.2.3)  b ut  with  "periodic  boundary  conditions"  defined  by 


u(o,t)  = u(L,t)  ; U^fcjt)  * ux(L,t) 


(1.3.24) 


(This  formulation  is  appropriate  to  the  longitudinal  vibrations  of  air  in 
a tube  closed  on  itself  (Fig.  1.4),  assuming  that  curvature  effects  are 
ignorable. ) 

For  the  exponential  function  (1.3.3)  it  is  readily  establis'.ed  that 


a = 2 (e-2)  , a = 2(e-l)  / (1  + 4n2ir2) 
o n 


(1.3.25) 


b = - (e-1)  4nir  / (l+4n2ir2) 


The  representation  obtained  by  summing  up  to  n = 4 is  shown  in  Fig.  1.3(d); 
since  the  extended  function  now  has  discontinuities  at  both  0 and  L,  the 
Gibbs  phenomenon  occurs  at  both  these  points. 

Sometimes  (1.3.21)  is  written  in  "polar  form",  namely. 


r(x)  = a /2  + Z r cos(2nnx/L  - <J>  ) 


n=l 


n 


(1.3.26) 


where 


n 


(a2  + b2)**  and  tan  $ = b /a 


n 


n n n 


(1.3.27) 


f 
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Note  that  here  orthogonality  involves  an  integral  of  the  form 


* 

/u  u dx 
m n 


(1.3.36) 


where  u and  u are  the  mth  and  nth  functions  of  the  set  and  the  asterisk 
in  n 

Note  that  the  derivation  of  the 


denotes  a complex  conjugate  value, 
coefficients  using  eqn  (1.3.35)  is  valid  when  g(x)  is  a complex  function. 
The  case  of  real  g(x)  is  characterised  by  the  relation 


(1.3.37) 


-n 


n 


This  follows  immediately  from  (1.3.34).  Parseval's  theorem  is  here 
established  by  multiplying  each  side  of  eqn  (1.3.33)  by  its  complex 
conjugate  and  integrating  over  (0,L).  The  result  is 


g g * dx  = (1/L) 


c c" 


n n 


(1.3.38) 


When  applied  to  the  exponential  function  (1.3.3),  the 


coefficients  c turn  out  to  be  given  by 
n 


e-2,  c = (e-1)  / (1-2-ir  in) 


(1.3.39) 


Sunning  the  complex  series  (1.3.33)  from  n = -N  to  n=N  gives  exactly  the 
same  result  as  summing  the  trigonometric  series  from  n=o  to  n=N.  It  must 
be  borne  in  mind  that  a pair  of  complex  terms  is  needed  to  give  one  of  the 
terms  in  the  series  (1.3.26);  this  is  sometimes  described  by  saying  that 
"a  positive  and  a negative  frequency  are  needed  in  the  complex  represent- 
ation to  give  a single  frequency  in  the  real  representation. 


1.4  The  Fourier  Integrals 


1.4.1  Fourier  Complex  Exponential  Integral 


All  of  the  Fourier  series  have  a periodic  character  and, 
as  already  seen,  throw  up  repetitions  of  the  original  function  defined  over 
the  basic  interval  (0,L).  However,  in  many  problems  where  the  independent 
variable  x ranges  from  — to  » there  is  a requirement  to  represent  a function 
which  is  non-periodic.  (e.g.  Referring  again  to  Fig.  1.3(a)  it  might  be 
required  to  represent  the  function  precisely  as  it  is  shown  there;  in 


* 


/ 


particular,  it  might  be  required  that  the  representation  returns  a zero 
value  of  the  function  for  all  x < o and  all  x > L.)  None  of  the  series 
can  do  this;  however,  it  earn  be  achieved  by  a Fourier  integral  as 
indicated  below. 


Consider  a function  g(x)  which  is  non-periodic  and  is 
essentially  negligible  outside  some  range  but  is  otherwise  more  or  less 
arbitrary.  Suppose  this  function  is  represented  by  a complex  Fourier 
series  of  period  L where  L more  than  spans  the  range  in  which  the  function 
is  non-zero  (Fig.  1.5a).  Here  it  is  convenient  to  take  the  range  as 
(-L/2,  L/2) , rather  than  (0,L)  as  previously,  and  it  is  easily  established 
that  the  representation  now  is  given  by 


c exp  (2ir  inx/L) 


with 


g(x)  exp  (-2tt inx/L)  dx 


This  representation,  of  course,  implies  the  periodic  repetitions  of  Fig. 1.5 (a) 
Now,  if  a larger  value  of  L is  chosen,  whilst  the  function  remains  unaltered, 
these  repetitions  will  be  spaced  further  apart  as  in  Fig. 1.5(b).  It  might 
be  hoped  that  in  the  limit  of  L approaching  infinity  the  repetitions  will 
be  lost  altogether;  this,  in  fact,  occurs  and,  concomitantly , the  sum  (1.4.1) 
becomes  an  integral.  This  can  be  seen  as  follows.  Think  of  L as  very  large 
but  (temporarily)  fixed.  If  a new  variable,  f,  defined  by 


is  introduced  then,  as  n ranges  over  all  integral  values  from  -°°  to  °°,  f will 
pass  through  a sequence  of  closely  spaced  values  from  to  00  . If  L is  large 

enough  f will  approximate  a continuous  variable,  and  the  difference  between 
successive  values  of  f can  be  written  as  the  differential,  df 


On  substituting  from  (1.4.3)  and  (1.4.4)  into  (1.4.2),  and  letting  L approach 
«,  it  follows  that 


! 


I 
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G(f)  - / g(x)  exp  (-2  irifx)  dx 


(1.4.6) 


Further,  on  substituting  from  (1.4.3)  and  (1.4.5)  into  (1.4.1),  and  on 
replacing  the  sum  of  infinitesimals  by  am  integral,  it  is  found  that 


g(x)  - / G (f ) exp  (2irifx)  df 


(1.4.7) 


The  two  results  (1.4.7)  and  (1.4.6)  represent  the  complex 
exponential  form  of  Fourier’s  integral  theorem;  they  are  the  respective 
amalogues  of  eqns  (1.3.33)  and  (1.3.34)  for  a non-periodic  function. 
Commonly,  G(f)  as  defined  by  (1.4.6)  is  said  to  be  the  "Fourier  transform" 
of  g(x)  and  the  relation  (1.4.7)  is  referred  to  as  the  "inversion  formula". 
(Instead  of  considering  the  Fourier  integral  as  applying  to  a non-periodic 
function,  sometimes  it  is  helpful  to  think  of  it  as  applying  to  a function 
with  an  infinite  period.) 


1.4.2 


Other  Fourier  Integrals 


The  Fourier  integral  just  discussed  was  derived  from,  and 
is  closely  analogous  to,  the  complex  form  of  the  Fourier  series.  In  a 
parallel  way  it  is  possible  to  derive  sine,  cosine  and  trigonometric  integrals 
which  are  analogous  to  the  other  series  forms.  For  details,  see  the 
references  already  cited.  In  all  that  follows,  however,  attention  will  be 
restricted  to  the  complex  form. 

1.5  Properties  and  Simple  Examples  of  Fourier  Transforms 


1.5.1 


Basic  Properties 


Here  some  of  the  basic  properties  of  Fourier  transforms, 
which  are  useful  in  manipulating  with  them  are  listed;  these  can  all  be 
readily  proven  from  the  definitions  (1.4.6)  and  (1.4.7).  Throughout,  the 
notation  of  "g"  for  the  function  and  "G"  for  its  transform  will  be  adopted. 


(a) 


Linearity  : If 


g = ci  gx+  °2  g2 


where  c^  and  are  constants,  then 


G " C1  °1  + C2  G2 


(1.5.1) 


(1.5.2) 


1 


Shift  Theorem  : If  gx  and  g2  are  such  that 


g2(x)  = gx  (x-a) 


where  a is  a constant, then 


G2(f)  - exp  (-2irifa)  G^f) 


Transform  of  Derivative  : If 


(1.5.3) 


(1.5.4) 


h(x)  = dg/  dx 


H = 2irifG 


(1.5.5) 


(1.5.6) 


(Here  it  is  assumed  that  g(x)  vanishes  "sufficiently  rapidly"  at 
infinity. ) 

Transform  of  Real  Function  : The  Fourier  transform  formulae  apply 
whether  g is  a real  or  complex  function.  If  g is  real,  then 


G (-f)  = G*  (f ) 


(1.5.7) 


Fourier  Transform  of  Rectangular  Pulse;  The  Sine  Function 


Consider  the  function  defined  by 


g(x)  = 1 


x|  <L/2 


x >L/2 


(1.5.8) 


This  function,  which  is  shown  in  Fig.  1.6,  is  widely  used  in  Fourier  analysis 
Direct  integration  readily  yields 


G (f ) = {sin(irf  L)  } /(iff) 


L sinc(irfL) 


Here  the  functional  notation 


sine  y = sin  y/y 


(1.5.9) 


(1.5.10) 


has  been  introduced.  (Sometimes  used  alternatives  to  "sine"  are  "sync" 
and  "dif"j  the  last  name  was  presumably  chosen  because  the  function  has  a 
fundamental  role  in  diffraction  problems.) 
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A sketch  of  the  function  sine  y is  shown  in  Fig.  1.7. 

The  function  has  a maximum  value  of  unity  at  y = 0 and  is  zero  at  y * i n » 
for  n = 1,2  •••.  Its  turning  points,  apart  from  the  origin,  are  approximately 
at  y = ± (2n  + l)x/2.  More  accurately,  the  turning  point  in  the  vicinity  of 
y = 3x/2(  =4.71)  is  at  y = 4.49;  the  other  turning  points  are  nearer  the 
approximate  values  quoted  above.  The  amplitude  of  the  minimum  at  y = 4.49 
is  0.22,  and  the  amplitudes  at  the  other  turning  points  essentially  are 
decreasing  as  1/y.  That  part  of  the  function  between  -it  and  iris  referred 
to  as  the  "main  lobe". 

Reverting  to  the  transform  (1.5.9),  it,  of  course,  also  is 
of  the  form  shown  in  Fig. 1.7;  here,  the  main  lobe  extends  from  f = -1/L  to 

f * 1/L.  Note  that  decreasing  L gives  a sharper  (or  more  concentrated) 
original  function  but  gives  a more  diffuse  (or  less  concentrated)  transform,  and 
conversely.  This  is  a general  feature  of  Fourier  transforms;  in  fact,  there 
exists  an  "uncertainty  principle"  (exactly  as  in  Quantum  Mechanics)  relating 
the  "width"  of  the  function  and  that  of  its  transform. 

Finally,  if  the  pulse  under  consideration  extends  from  0 to 
L,  rather  than  from  -L/2  to  L/2,  the  "shift  theorem"  4s  given  by  eqn  (1.5.4) 
can  be  used  to  write  down  the  transform.  The  result  is  now 

G(f)  = L exp  (-it  ifL)  sinc(irfL)  (1.5.11) 

1.5.3  Fourier  Transform  of  Triangular  Pulse 

Now  consider  the  "triangular"  function  defined  by 

g (x)  = (1  + 2x/L)  -L/2  < x < 0 

= (1  - 2x/L)  0 < x < L/2  (1.5.12) 

=0  |KI  > L/2 

This  function  is  shown  in  Fig.  1.8;  it  is  also  of  common  occurence  in  Fourier 
analysis.  Only  elementary  integration  is  needed  to  establish  that  the 
transform  is 

G (f ) = (2/L)  {sin2  (irfL/2)  } /(ir2f2) . (1.5.13) 

It  can  be  seen  that  the  behaviour  of  the  transform  is  essentially  that  of 
the  function 

sin2  y / y2  (1.5.14) 


i 
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This  last  is  shown  in  Fig.  1.9.  It  is  of  a generally  similar  character  to 
the  sine  function,  the  most  important  differences  being,  firstly,  that  die 
present  function  is  always  positive  and,  secondly,  there  is  a much  greater 
concentration  of  the  function  in  the  "main  lobe". 

1.6  The  Dirichlet  Kernel 

The  arguments  given  to  date  for  the  development  of  both  Fourier  series 
and  integrals  are  more  or  less  in  accord  with  what  took  place  historically. 
However,  these  arguments  have  been  plausible,  rather  than  rigorous.  The 
rigorous  proof  of  the  Fourier  expansions  was  given  by  Dirichlet  in  1829  and 
an  outline  of  this  proof  is  given  below  for  the  Fourier  integral.  (The  proof 
for  the  series  is  very  similar) . 

The  Fourier  integral  theorem  is  embodied  in  the  dual  relations 

g(x)  =/°°  G (f ) exp (2irifx)  df  (1.6.1) 

G (f ) */°°  g(x)  exp(-2irifx)  dx  (1.6.2) 

If  these  relations  are  valid  then  on  substituting  from  (1.6.2)  into  (1.6.1) 
an  identity  should  result.  It  is  convenient  to  commence  by  taking  finite 
integration  terminals  in  (1.6.1),  say  -N  and  N where  N is  large;  the 
associated  function  on  the  left  hand  side  of  (1.6.1)  will  now  be  written  as 
gN(x)  i.e., 

g„(x)  = /*!  G (f ) exp (2irifx)  dx  (1.6.3) 

Thus,  it  is  required  to  show  that  as  N approaches  infinity,  gN  approaches 
g.  Substituting  for  G from  (1.6.2)  - where  u is  now  written  for  the 
integration  variable  - into  (1.6.3)  and  reversing  the  order  of  integration 
gives 

g„(x)  = / °°g(u){  / *?  exp  2irif (x-u)  df}  du  (1.6.4) 

N -00  ~N 

i.e., 

gN(x)  = dn<x-u>  du  (1.6.5) 


where 
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D (x-u)  « sin  { 2IIN  (x-u)  }/  n(x-u) 
N 


2N  sinc{  2I1N (x-u)} 


(1.6.6) 


(1.6.7) 


If  (x)  is  to  approach  g(x),  then  the  limiting  value  D(x-u)  of  the  kernel 
function  DN(x-u)  must  have  the  infinitely  strong  "focussing  power"  implied 
by  the  equation 


g(x)  = / g(u)D(x-u)  du 

—00 


(1.6.8) 


The  limit  function  D(x-u)  is  termed  the  "Dirichlet  kernel".  (From  (1.6.7) 
it  is  clear  that  the  "main  lobe"  of  the  sine  function  in  this  kernel  is  centred 
around  *=u;  also  that  the  amplitude  of  the  kernel  becomes  infinite  in  the  limit.) 


Introducing  a new  variable  y defined  by 


y * 2IIN  (x-u) 

it  is  easily  seen  that  (1.6.5)  becomes 

g (x)  = (1/TI)  / °°g(x  - y / (2IIN)  ) sine  y dy 
N — 00 


(1.6.9) 


(1.6.10) 


If  N now  is  let  approach  infinity,  and  if  use  be  made  of  the  standard  integral, 

/ °°sinc  y dy  = II  (1.6.11) 

—•CO 

it  follows  that  gN(x)  approaches  g(x)  as  required. 

The  above  result  is  valid  as  long  as  g(x)  is  a continuous  function 
(and  becomes  suitably  small  at  infinity.)  If  g(x)  has  a point  of  discontinuity, 
then  a more  rigorous  analysis  shows  that,  at  such  a point,  gN(x)  approaches  the 
mean  of  the  values  on  either  side  of  the  discontinuity. 

Before  passing  on  it  might  be  observed  that  the  Dirichlet  kernel 
obviously  has  an  intimate  connection  with  the  sine  function  which,  in  turn, 
appeared  as  the  transform  of  a rectangular  pulse.  Another  kernel  function 
sometimes  encountered  is  the  "Fejer  kernel";  this  has  even  stronger  focussing 
powers  than  the  Dirichlet  kernel  and  is  closely  connected  with  the  transform 
of  a triangular  pulse. 
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1.7  Delta  Functions;  Parseval1 s Theorem  for  Fourier  Integral 


1.7.1  General 

Consider  again  the  Dirichlet  kernel  D(x-u)  just  discussed. 
This  was  defined  as  the  limit  for  N approaching  infinity  of 

D„(x-u)  = 2N  sine  {2  JlN(x-u)}  (1.7.1) 

N 

The  kernel  has  the  following  properties  : 

f °°D (x-u)  du  = 1 (1.7.2) 

*00 

/_* g(u)  D(x-u)  du  = g(x)  (1.7.3) 

The  first  of  these  properties  follows  from  the  fact  that  (1.7.2)  holds 
when  D,  replaces  D,  for  all  N.  (This  can  be  proved  by  direct  integration.) 
Hence  it  is  also  true  for  the  limit  function.  The  second  property  has 
already  been  demonstrated  in  the  previous  section. 


The  relation  (1.7.3)  shows  that  D can  be  regarded  as  a 
"weighting  function"  which  gives  a zero  weight  to  a function  g(u)  at  all 
points  of  an  infinite  interval  save  at  the  single  point  u=x,  where  it 
assigns  a unit  weight.  Ordinary  weighting  functions  will  not  possess  this 
property.  Functions  which  do  are  variously  termed  "delta  functions", 
"distributions",  "impulse  functions"  or  "generalised  functions".  Detailed 
studies  of  their  properties  can  be  found,  for  instance,  in  refs.  11  and  12. 


Actually  the  Dirichlet  kernel,  although  the  commonest  type 
of  delta  function  occurring  in  Fourier  analysis,  is  somewhat  more  complicated 
than  other  delta  functions.  The  simplest  such  function  is  that  defined  by 


6 (x-u)  = lim  e-K)  1/e  for  x < u < X +£ 
= o all  other  u 


(1.7.4) 


This  function  is  shown  in  Fig.  1.10  for  (a)  non-zero  s and  (b)  zero  e.  One 
physical  interpretation  of  it  is  as  a band  of  uniform  (line)  pressure  which 
in  the  limit  becomes  a concentrated  force.  It  is  immediately  established 
that  a relation  analogous  to  (1.7.2)  holds.  To  prove  the  result  analogous  to 


I 
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(1.7.3)  suppose  that  g(u)  is  expanded  in  a Taylor  series  about  the  point  x. 
Then,  writing  v ■ u-x,  it  follows  that 

I "g(u)  6 (x-u)  du  = lim  e-HD  /*  (1/e){  g(x)  + vgX(x)  + •••Jdv 


lira  fc-K)  {g(x)  +£g1(x)/2  + •••} 


i . e . , 


f g (u)  6 (x-u)du  = g(x) 


(1.7.5) 


However,  as  already  indicated,  in  most  of  the  following 
the  delta  function  used  will  be  that  associated  with  the  Dirichlet  kernel. 
This  can  be  expressed  in  the  form  (using  the  "<S"  notation) 


/ °°  exp  {2iri  (f  -f)  x}  dx  = Mf.-f) 

— OO  X X 


(1.7.6) 


Also,  it  might  be  noted  that  since  the  sine  function  is  an  even  one,  it 
follows  that  here  at  any  rate 


fitfj-f)  = Mf-f^  (1.7.7) 


1.7.2  Fourier  Transforms  involving  Delta  Functions 

Conventional  Fourier  transforms  can  only  be  applied  to 
functions  g(x)  for  which  some  boundedness  condition  such  as 


I_K  l ^ (x) | dx  < M 


(1.7.8) 


where  M is  finite,  applies.  However,  using  delta  functions,  they  can  be 
fruitfully  applied  to  common  functions  for  which  (1.7.8)  does  not  hold. 


(a)  Transform  of  constant  : Consider  the  case  where 


i 


g(x) 


all  x 


(1.7.9) 


(This  can  be  regarded  as  the  limit  of  the  rectangular  pulse 
discussed  in  Section  1.5.2  when  the  length  L of  the  pulse  approaches 
infinity. ) 
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The  transform  of  (1.7.9)  is 


G(f)  = / exp  (-27rifx)  dx  = 6(f) 


(1.7.10) 


i.e.,  the  transform  is  a delta  function  located  at  the  origin  of 
f.  The  function  and  its  transform  are  shown  in  Fig.  1.11  and  it 
can  be  seen  that  this  represents  an  extreme  case  of  the 
"uncertainty  principle"  already  referred  to  in  Section  1.5.2. 

Transform  of  periodic  function  : As  a second  example/  consider 


g(x)  = exp  (2irinx/L) 


(1.7.11) 


where  n is  an  integer;  this,  of  course,  is  the  typical  term  in 
the  complex  Fourier  series  expansion  of  a function  of  period  L. 
The  transform  is  easily  found  to  be 


G (f ) = 6(  f - n/L  ) 


(1.7.12) 


Since  an  arbitrary  periodic  function  can  be  represented  as  an 
infinite  series  of  terms  such  as  (1.7.11),  with  n ranging  through 
all  positive  and  negative  integral  values,  its  transform  will 
consist  of  an  infinite  number  of  delta  functions  located  at  the 
points  f = n/L. 


If 


g(x)  = cos(2irnx/L) 


(1.7.13) 


then,  immediately, 


G (f ) = (1/2)  {6 (f-n/L)  + 6 (f+n/L)} 


(1.7.14) 


This  function  and  its  transform  are  shown  in  Fig.  1.12  for  the  case 
n = 2. 


As  a final  example  of  the  use  of  delta  functions,  Parseval’s 
theorem  for  the  Fourier  integral  will  be  proved. 


As  always. 


g(x)  = / G (f ) exp(2irifx)  df 


(1.7.15) 


Taking  the  conjugate  complex  of  each  side,  and  writing  fj  for  the 
integration  variable,  gives 


g*  (x)  = / G* (f i ) exp(-2irif jx)  dfj 


(1.7.16) 


Hence, 


/"g  g*  dx  = /*G(f)  {/“g*  (fj){  /“exp  2iri  (f-fi)  x dx  >df df 

(1.7.17) 

The  innermost  integral  is  simply  6(f-fi)  and  so  the  middle  integral  returns 
only  G*  (f ) ; the  final  result  is 


/"g  g*  djc  = /„ G G* 


This  is  Parseval's  theorem. 


1.8  Convolutions 


(1.7.18) 


An  inportant  tool  in  all  Fourier  analysis  is  the  convolution  integral. 
Its  introduction  can  be  motivated  in  a variety  of  ways.  One  way  is  via  the 
following  question  : Suppose  that  the  Fourier  transform  of  seme  function 
g(x)  is  sought  and  the  direct  calculation  of  the  integral  is  difficult  (which 
is  more  often  the  case  than  might  be  thought  frem  its  rather  simple  form) , but 
suppose  that  the  function  is  recognized  as  the  product  of  two  other  functions 


whose  transforms  are  known  i.e.. 


g(x)  = gi (x)  g2<*> 


(1.8.1) 


where  Gj (f ) and  G2<f)  are  known.  Can  the  transform,  G,  of  g be  expressed 
in  terms  of  Gi  and  G2  ? 

From  the  definition. 


G (f ) = / gi  (x)  g2(x)  exp  (-2irifx)  dx 


(1.8.2) 


Also, 


gi  (x)  = Gi  (f  1 ) exp  (2TTiflX)  dfx 


g2(x)  = / G2(f2)  exp  (2irif2X)  dfz 


(1.8.3) 


(1.8.4) 


On  substituting  from  (1.8.3)  and  (1.8.4)  into  (1.8.2),  and  reversing  the 
order  of  integration  so  that  the  x integration  is  done  first,  it  is  found 


G(f ) = /^Gi  (f  1 ) { /“g2  (f  2 ) t /”exp(2iri  (fx  + f2  -f)  x)  dx}df2}dfx 

(1.8.5) 


The  innermost  integral  can  be  written  as 


6 ( f2  (f-fi)  ) 

and  this,  when  applied  to  the  middle  integral,  gives  simply 


G2  ( f ~f 1 ) 

Hence  the  final  result  is 


G (f)  = fjZltfl)  G2 (f~f 1 ) df j 


(1.8.6) 


An  integral  of  this  form  is  termed  a "convolution"  of  the  functions  Gx  and 
G2  • Sometimes  it  is  denoted  by  a special  symbol  e.g.. 


G — Gj*  G2 


(1.8.7) 


(It  is  easily  established  that  it  does  not  matter  which  of  Gj  and  G2  has 
the  argument  "f j"  and  which  " f-f j" . ) 
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A quite  analogous  proof  can  be  developed  to  show  that,  if 
G(f)  - G^f)  G2(f)  (1.8.8) 


then 

g(jc)  = /"gi(*i)  g2  (x-xi)  dxi  (1.8.9) 

As  a simple  example  of  these  last  relations  the  "shift  theorem"  of 
eqn  (1.5.4)  will  be  proved  in  reverse.  Suppose 

G(f)  = exp  (-2irifa)  H (f)  (1.8.10) 

where  h(x)  is  known  and  g(x)  is  sought.  This  is  a product  of  the  form  of 
(1.8.8).  The  inverse  of  the  exponential  term  is 

/ooexp(-2nifa)  exp(2irifx)  d f = 6(x-a)  (1.8.11) 

Hence,  by  the  convolution  theorem  (1.8.9), 

g(x)  = /“  6 (xj-a)  h ( x-Xl)  dxj 

i.e. , 

g (x)  = h(x-a) 
as  required. 

1.9  Fourier  Expansions  in  Two  or  More  Variables 

All  the  discussion  to  date  has  related  to  Fourier  expansions  of  a 
function  of  one  variable.  However  the  concepts  can  be  readily  extended 
to  functions  of  two  or  more  variables.  This  will  be  demonstrated  for  the 
Fourier  integral  but  a similar  approach  can  be  used  for  the  series. 

t 

Consider  a function  g(xj,  X2)  and,  temporarily  thinking  of  X2  as 
fixed,  take  its  transform  with  respect  to  xj  to  get 

61  (fi,  x2)  - /”g(xlf  x2)  exp  (-2-nifjXi)  dx! 


(1.8.12) 

(1.8.13) 


(1.9.1) 


with  the  inverse  relation 


g(x1#x2)  = Gi(fi,  x2 ) exp  (2iri  fjx^  dfi 


Now  take  the  transform  of  Gj  with  respect  to  x2 


to  get 


(1.9.2) 


G(fi,  f 2)  - Gi(fi,  x2)  exp  (-2irif2x2)  dx2 


(1.9.3) 


with  the  inverse  relation 


Gl  (fi,x2)  = exp  (2xif2x2)  df2 


(1.9.4) 


Substituting  from  (1.9.1)  into  (1.9.3)  gives 

G(f  1 *f2)  * /“  /”g(xi,  x2)  exp  {-2iri(f1x1  + f2x2)}  dXldx2 

(1.9.5) 

and  from  (1.9.4)  into  (1.9.2)  gives 

g(x1#  x2 ) = /*  /°*G (f  1 » f2)  exp  {2iri(f1x1+  f2  x2)  dfjdf2 

(1.9.6) 

The  two  relations  (1.9.5)  and  (1.9.6)  represent  the  Fourier  integral  theorem 
for  a function  of  two  variables.  The  extensions  to  functions  of  more  than 
two  variables  cure  obvious. 
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APPENDIX 

End  Discontinuity  Effects  and  Hannin 


In  the  earlier  "physical"  derivation  of  the  different  Fourier 
series  it  transpired  that  each  form  had  its  own  particular  boundary 
conditions  associated  with  the  end-points  0 and  L of  the  fundamental 
interval.  Again  writing  g(x)  for  the  function  to  be  represented,  these 
boundary  conditions  were  as  follows  : 


sine  series  (period  2L) 


cosine  series  (period  2L) 


, g(o)  = g(L)  = 0 


, g1 (o)  = g1 (L)  - 0 


trigonometric  series  (period  L)  , g(o)  = g(L)  and  g1  (o)  = g1  (L)  (A3) 

These  conditions  are  not  necessary  for  the  validity  of  the  series 
expansions;  however,  from  the  practical  viewpoint,  the  series  will  generally 
converge  much  more  rapidly  when  the  conditions  are  fulfilled. 

In  some  situations  the  following  device  is  useful.  Suppose  g(x) 
is  the  function  which  it  is  wished  to  represent,  but  suppose  that  it  does 
not  satisfy  the  "physical  boundary  conditions".  Let  h(x)  be  another 
function,  which  has  the  same  periodicity  as  the  expansion  to  be  used,  and 
which  is  chosen  such  that  the  "modified  function" 


m(x)  = g (x)  h(x)  (A4) 

does  have  the  physical  boundary  conditions.  Then  m(x)  can  be  expanded  in 
the  same  type  of  series  as  was  planned  for  g(x)  and  this  series  should  have 
good  convergence  properties. 

For  example,  suppose  a sine  series  is  to  be  used  and  that  g(o)  =0 
as  required  (as  would  be  the  case  for  an  odd  function)  but  that  g(L)  f 0. 


Here  taking 


(1  + cos  irx/L)/2 


satisfies  the  requirement  that  h(L)  =0  and,  also,  h(x)  is  of  period  2L. 

This  function,  which  is  known  as  a Hann  function  (after  J.  von  Hann)  is  shown 

in  Fig.  Al.  If  the  Fourier  coefficients,  b , for  the  sine  series  for  g(x) 

n 

are  known,  the  coefficients  b^  for  the  modified  function  m(x)  can  be  written 
down.  From  the  definitions. 
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b » (2/L)  / g(x)  sin(ntrx/L)  dx 
n o 


(A6) 


b * (2/L)  / g(x)h(x)sin(nirx/L)  dx 
n o 


(A7) 


With  h (x)  as  given  by  eqn  (A5) , it  follows  that 


h (x) sin (nux/L)  = (1/4)  sin{  (n-l)irx/L}  + (1/2) sin(nirx/U 

+ (1/4)  sin  { (n+l)irx/L> 


(A8) 


Hence , 


bi/2  + b2/4 


(A9) 


b„  -b  ./4  + b /2  + b /4  , n*2, 3 

n n-l  n n+i 


As  an  example,  consider  the  Fourier  sine  series  expansion  for 
the  function 


g(x)  * exp  (x/L)  - 1 (A10) 

This  has  already  been  determined  in  Section  1.3.1.  The  values  of  the  first 
few  coefficien 
Table  1 below. 


few  coefficients,  b , and  the  corresponding  "harmed"  ones,  b*,  are  given  in 
n n 


n 

b 

n 

b1 

n 

1 

0.876 

0.305 

2 

-0.533 

0.041 

3 

0.356 

-0.023 

4 

-0.272 

0.007 

Table  1 : Fourier  Coefficients  for  Exponential 
Function  before  and  after  Hanning. 


i 


Clearly  the  harmed  coefficients  are  decreasing  much  more 
rapidly  so  that  taking  only  4 terms,  say,  in  the  series  for  the  hanned 
function  should  give  a much  better  representation  of  it,  than  does  the 
same  number  of  terms  in  the  series  for  the  original  function.  Thus,  if 
one  is  faced  with  the  problem  of  recovering  a function  given  its  Fourier 
sine  coefficients  the  following  procedure  may  be  useful  : 


(i)  From  the  original  coefficients,  determine  the  hanned 
coefficients  using  eqn  (A9) . 

(ii)  Recover  the  modified  (or  hanned)  function  by  summing 
the  series  in  the  hanned  coefficients. 

(iii)  Recover  the  original  function  from  the  hanned  one  using 
g = m/h. 

The  results  of  doing  this  for  the  exponential  function  (A10) 
are  shown  in  Table  2 below  (using  4 terms  in  the  series) . For  comparison 
the  exact  results  and  the  results  from  a direct  summation  of  the  series 
(again  with  4 terms)  cure  included. 


x/L 

g(x)  = exp  (x/L)  - 1 

Exact 

Direct  Sum 
(n=4) 

Hanning 

(n=4) 

0 

0 

0 

0 

0.1 

0.105 

-0.013 

0.109 

0.2 

0.221 

0.186 

0.222 

0.3 

0.350 

0.471 

0.345 

0.4 

0.492 

0.569 

0.490 

0.5 

0.649 

0.520 

\ 

0.656 

0.6 

0.822 

0.679 

0.829 

0.7 

1.014 

1.166 

0.991 

0.8 

1.226 

1.521 

1.186 

0.9 

1.460 

1.131 

1.800 

1.0 

1.718 

0 

- 

Table  2 : Hanning  with  Fourier  Sine  Series  for  Exponential 
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When  dealing  with  a Fourier  cosine  series  for  a function  which 
has  zero  slope  at  x ■ 0 but  not  at  x = L exactly  the  same  device  can  be 
used.  It  can  be  readily  established  that  the  modified  function  (A4)  with 
h(x)  again  given  by  (A5)  will  have  a zero  slope  at  both  end  points.  The 
hanned  coefficients,  a1,  are  given  in  terms  of  the  original  ones,  a , by 


a /2  + 
o 


a /2 


a /4  + a /2  + a ./4 
n-1  n n+1 


(All) 


Finally,  when  dealing  with  the  general  trigonometric  series  for 
a function  of  period  L,  but  which  does  not  satisfy  the  periodic  boundary 
conditions,  use  of  the  Hann  function 


h(x)  * (1-cos  2irx/L)/2 


(A12) 


leads  to  a modified  function,  m(x),  which  will  have  periodic  boundary 

conditions.  (In  fact,  m(x)  will  be  zero  at  x*0  and  x=L,  as  will  be  its 

derivative).  This  Hann  function  is  shown  in  Fig.A2.  By  a procedure 

analogous  to  that  used  in  establishing  eqn  (A8)  it  can  be  shown  that  the 

hanned  coefficients  a1,  b1  are  related  to  the  original  ones  a , b as 

n n n n 

follows  : 


a /2  - a./2 
o 1 


n 


"Vi/4  + V2  - Vi/4 


(A13) 


bj  = bj/2  - b2/4 


- b _/4  + b/2  - b , / 4 
n-1  n n+i 


The  procedure  for  using  hanning  here  is  completely  analogous  to 
that  already  described  for  the  sine  series. 

The  above  is  not  the  only  method  for  dealing  with  the  problem 
end  effects;  another  way  is  to  subtract  off  a "linear  trend". 


(d)  Trigonometric  series  representation 


O Numerical  evaluation  of  series,  summing  to  n = 4 

Function  represented  by  (infinite)  Fourier  series 

outside  basic  interval 


FIG.  1.3  VARIOUS  FOURIER  SERIES  REPRESENTATIONS 
OF  THE  ONE  FUNCTION:  g(x)  * EXP(x/L)  - 1 


FIG.  1.4  TUBE  OF  AIR  ILLUSTRATING  PERIODIC  BOUNDARY 
CONDITIONS  (CURVATURE  EFFECTS  IGNORED  SO 
ONE-DIMENSIONAL  WAVE  EQUATION  APPLIES) 


-5L/2  -3L/2  -L/2  0 L/2  3L/2 

(a)  Periodic  function,  basic  interval  L 


(b)  Same  periodic  function,  basic  interval  L1  = 2L 


FIG.  1.5  EFFECT  OF  INCREASING  BASIC  INTERVAL 
ON  REPETITIONS  OF  PERIODIC  FUNCTION 
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CHAPTER  2 : THE  DISCRETE  FOURIER  TRANSFORM 

2.1  General 

The  analytical  processes  associated  with  the  classical  Fourier 
series  and  integrals  can  only  be  conveniently  performed  when  the  functions 
considered  are  of  a relatively  simple  form.  In  other  cases,  and,  in 
particular,  for  functions  that  are  only  specified  numerically  (as  is  the 
general  situation  when  the  functions  are  experimental  data)  it  is 
necessary  to  adopt  numerical  Fourier  methods.  These  have  a history  as  long  as 
the  analytical  methods,  going  back  to  Lagrange  who,  around  1760,  and  as  a 
result  of  studies  of  the  vibration  of  a string  with  particles  attached, 
constructed  a theory  of  interpolation  using  trigonometric  functions  which 
is  quite  analogous  to  the  more  familiar  theory  based  on  polynomials. 

Lagrange's  results  form  the  basis  of  what  today  are  called  "discrete  Fourier 
transforms"  (DFTs) . 

Sporadic  developments  in  this  area  occurred  over  the  next  two 
centuries  or  so.  For  example,  one  problem  to  which  Lagrange's  results  were 
applied  was  in  the  search  for  hidden  periodicities  in  experimental  data;  an 
important  tool  here  was  the  "periodogram"  developed  by  Schuster  around  1900. 
Also,  Runge  around  the  same  time  developed  computational  procedures  for  use 
in  trigonometric  interpolation  which  embodied  the  key  features  of  the  "fast 
Fourier  transform"  (FFT)  used  today.  Some  account  of  this  early  work  is 
given  in  ref.l.  Danielson  and  Lanczos  in  1942  treated  several  important 
aspects  of  numerical  Fourier  analysis;  this  work  is  elaborated  on  by  Lanczos 
in  reference  2 and  contains  the  seeds  of  many  present  day  developments.  In 
the  post  World  War  2 era  the  whole  subject  has  been  greatly  advanced, 
especially  by  workers  in  the  field  of  communication  engineering.  (See,  for 
exanqple,  ref. 3)  This  culminated  in  1965,  with  Cooley  and  Tukey  developing 
the  definitive  form  of  the  fast  Fourier  transform  which  provided  a dramatic 
reduction  in  the  computing  effort  previously  required. 

Nevertheless,  amongst  the  spate  of  books  on  general  numerical 
analysis  which  appeared  following  the  advent  of  the  digital  computer,  very 
few  made  more  than  a fragmentary  reference  to  numerical  Fourier  methods; 
an  exception  was  that  by  Hamming  (4).  (See,  also,  ref. 5).  It  is  only  over 
the  past  two  or  three  years  that  these  methods  have  been  treated  in  detail 
in  the  t*xt  book  literature  (e.g.  ref. 6). 


r 
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2.2  Physical  Motivation  for  Discrete  Fourier  Transform 

Just  as  the  vibrations  of  a continuous  (heavy)  string  provided 
the  physical  basis  for  the  development  of  Fourier  series,  so  do  the 
vibrations  of  a (light)  string  with  discrete  masses  attached  provide 
a physical  basis  for  the  development  of  the  DFT.  As  already  remarked 
this  problem  was  solved  by  Lagrange,  and  details  can  be  found  in  many 
texts  e.g.  that  by  Hildebrand  (7) . Only  the  outlines  of  the  solution 
will  be  given  here. 

Consider  the  case  of  N + 1 equal  particles  uniformly  spaced  on 
a string  of  length  L (Fig.  2.1).  The  particles  are  numbered  "0"  to  "N" 
as  shown  and  clearly  the  base  co-ordinate  of  particle  k is  given  by 


= k L/N 


(2.2.1) 


Analogously  to  the  continuous  string,  first  consider  the  situation  when 

the  string  .has  fixed  ends  (so  that  the  two  end  particles  do  not  participate 

in  the  motion)  and  commences  vibrating  with  zero  initial  velocity  from  some 

deformed  configuration  (Fig.  2.2).  The  initial  position  of  the  m th 

particle  is  denoted  by  g . It  is  convenient  to  use  vector  notation  in  the 

xn 

following  and  to  write,  for  example, 


(2.2.2) 


Similarly  using  u to  denote  the  set  of  particle  displacements  at  time  t 
it  can  be  shown  that  the  solution  of  the  (N-l)  ordinary  differential 
equations  which  govern  this  problem  is  given  by  , 


U(t)  = l B e cos  p t 
, n n n 

n=l 


(2.2.3) 


where  the  Br  are  integration  constants,  the  p^  cure  the  natural  frequency 
parameters  (whose  precise  values  are  obtained  but  are  not  of  concern  here) , 
and  the  e^  are  the  normal  modes  given  by 


* 
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(2.2.4) 


e 

n 


sin  (rni/N) 
sin  (2rot/N) 


sin  ((N-l)nir/N)J 


(Bearing  in  mind  eqn  (2.2.1)  it  can  be  seen  that  (2.2.4)  represents  an 
approximation  to  the  n th  mode  for  the  continuous  string  which  was 
sin(nirx/L)  .)  The  expression  (2.2.3)  satisfies  the  condition  of  zero 
initial  velocity  whilst  the  requirement  that  u(o)  = g implies 


9 


N-l 


£ 

n=l 


B 

n 


e 

n 


(2.2.5) 


Since  equality  of  a vector  relation  implies  equality  of  corresponding 
components  eqn  (2.2.5)  can  just  as  well  be  written 

N-l 

g = £ B sin(  mnir/N)  (2.2.6) 

m , n 
n=l 

Since  the  g's  are  arbitrary,  eqn  (2.2.6)  shows  that  an  arbitrary  function 
defined  on  a discrete  set  of  uniformly  spaced  points  can  be  represented  by 
a finite  series  of  sines. 

As  in  the  continuous  case,  the  normal  modes  are  orthogonal  but 

here  orthogonality  is  expressed  by  the  vanishing  of  the  scalar  product  of 

e and  e , say;  i.e.,  using  the  bracket  notation  for  the  scalar  product, 
n k 

^ N-l 

(e  ,e,  ) = £ sin  (mnir/N)  sin(mkir/N)  =0,  n / k (2.2.7) 

n k ___  • 


For  an  analytical  proof  of  this  result,  see  ref.  7,  where  it  is  also  shown 
that 


(e  ,e  ) = N/2 

n n 


By 

virtue  of  these  relations, 

> “►  . 

B 

= (g,  e ) / 

(e  , e ) 

n 

n 

n n 

N-l 

B 

= (2/N)  £ 

g sin(mnir/N) 

n 

m=l 

m 

(2.2.8) 

it  follows  that 

(2.2.9) 

(2.2.10) 


\ 


The  dual  relations  (2.2.6)  and  (2.2.10)  constitute  a discrete  Fourier  sine 
transform  (DFST)  pair.  (Just  as  a note  of  historical  interest  it  might  be 
remarked  that,  whilst  Lagrange  derived  the  above  results  correctly,  he  then 
made  an  error  when  he  proceeded  to  the  limit  of  a continuous  string;  this 
led  him  to  pronounce  against  Bernoulli's  work.) 

Finally,  on  taking  the  scalar  product  of  each  side  of  eqn  (2.2.5) 
with  itself  the  analogue  of  Parseval's  theorem  is  obtained  in  the  form 


N-l  N-l 

E g2  = (N/2)  S B* 
m=l  n=l 


(2.2.11) 


Continuing  on  a parallel  argument  to  that  adopted  for  the  uniform 
string,  consider  the  case  where  the  present  string  has  "free  ends".  Now, 
the  two  end  particles  will  participate  in  the  motion.  It  turns  out  that 
the  result  analogous  to  eqn  (2.2.6)  is 


A cos(mnir/N) 
n 


(2.2.12) 


where  m is  in  the  range  (0,N)  and  the  two  accents  on  the  summation 
indicate  that  the  first  and  last  terms  are  each  taken  with  a weight  factor 
of  1/2  e.g.  the  first  term  in  (2.2.12)  is  Ao/2.  Again,  there  exist 
orthogonality  relations  for  the  discrete  cosine  functions  analogous  to 
eqns  (2.2.7)  and  (2.2.8),  and  these  can  be  used  to  show  that 


A = (2/N)  E g cos  (mnir/N) 

n m 

m=o 


(2.2.13) 


Since  the  g's  are  arbitrary,  eqn  (2.2.12)  shows  that  an  arbitrary 
function  defined  on  a discrete  set  of  uniformly  spaced  points  can  be 
represented  by  a finite  series  of  cosines.  (The  reason  for  the  weight 
factors  of  1/2  on  the  first  and  last  terms  is,  broadly,  that  in  considering 
the  discrete  case  as  an  approximation  to  the  continuous  case,  the  first  and 
last  particles  are  associated  with  a sub-interval  that  is  only  half  the 
length  of  that  associated  with  an  interior  particle.) 

Parseval's  theorem  here  takes  the  form 


(N/2)  E 


(2.2.14) 


l 


I 


2.3 


The  Different  DFTs 


2.3.1  The  DFST 


The  previous  discussion  has  shown  that  a function,  g, 
defined  at  a discrete  set  of  uniformly  spaced  points  over  an  interval 
(0,L)  can  be  represented  as  a finite  sine  series  : 


N-l 

E B sin  (mnir/N) 
, n 


where 


B 

n 


N-l 

(2/N)  E g sin(mnif/N) 
, m 
m=l 


(2.3.1) 


(2.3.2) 


Also,  it  will  be  recalled, 

x = m L/N 
in 


(2.3.3) 


m 


= 9(xJ 

in 


(2.3.4) 


As  with  the  infinite  series,  it  is  useful  to  look  at  the  function  which 
is  represented  by  the  finite  series  at  points  outside  the  basic  interval 
(0,L).  From  (2.3.3),  negative  values  of  x correspond  to  negative 
values  of  m and,  since  the  sine  is  an  odd  function,  the  series  (2.3.1) 
returns  the  result 

g = -g  (2.3.5) 

— m m 


Also,  since 

sin  (mnir/N)  = sin  { (m+2N)nir/N  } (2.3.6) 

the  series  (2.3.1)  represents  a periodic  function  of  period  2N.  These 
results  are  completely  analogous  to  the  continuous  case.  (Note  that  2N 
points  span  a length  of  2L. ) 
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2.3.2.  The  DFCT 

Likewise,  the  point  function  g can  be  represented 
as  a finite  cosine  series  : 


n=o 


A cos  (mnir/N) 
n 


where 


A 

n 


11 

(2/N) 

m=o 


cos  (mnir/N) 


(2.3.7) 


(2.3.8) 


The  series  (2.3.7)  is  again  of  period  2N  but  now  returns  a 
result  of  the  form 


9 


-m 


m 


i.e.  the  series  represents  an  even  function. 


(2.3.9) 


2.3.3.  The  DFTT 

If  it  is  required  to  represent  a function  g on  the 
same  set  of  points  as  in  the  above  but  now  witn  the  requirement  that, 
outside  the  basic  interval  (0,L) , the  extended  function  is  defined  by 
the  requirement  that  it  has  a period  of  N (or  L)  it  is  necessary  to 
use  a finite  series  of  combined  sines  and  cosines.  This  may  be 
designated  the  discrete  Fourier  trigonometric  transform  (DFTT) . The 
argument  is  basically  parallel  to  the  infinite  series  case. 

The  periodicity  requirement  can  be  expressed  as 


gm 


gm+N 


(2.3.10) 


in  particular  note  that 


(2.3.11) 


i 


« 


I 


so  that  here  one  can  only  specify  N arbitrary  values  of  g.  (In  the 
DFCT  one  could  also  specify  gN  arbitrarily;  in  the  DFST , gQ  and  gN 
were  implicitly  assumed  to  be  zero.) 


As  in  the  continuous  function  case,  write  g as  the  sum 
of  an  even  and  an  odd  point  function  : 


E + 0 
m m 


(2.3.12) 


where 


E = (1/2)  (g  + g ) 
xn  m —in 


(2.3.13) 


0 = (1/2)  (g  - g) 

v m m — m 


(2.3.14) 


Both  E and  o will  be  of  period  N;  hence,  these  last  two  relations  can 
be  just  as  well  written  as 


E = (1/2)  (g  + g ) 
m m n— m 


°m  " <1/2)  <*m  * W 


(2.3.15) 

(2.3.16) 


For  the  convenience  it  will  be  assumed  that  N is  an  even  number.  (As 

will  be  seen  later,  in  practice  N is  generally  a power  of  2.)  In  Fig. 2. 3, 

E and  0 are  shown  for  the  point  function  g derived  from  the 
mm  m 

continuous  function 


■g  (x)  = exp  (x/L)  -1 


(2/3. 17) 


talking  N=10.  It  will  be  noted  that  E is  symmetric  about  the  midpoint 

m 

m = N/2  and  that  0 is  skew-symmetric  about  that  point.  These  results 
m 

are,  of  course,  implied  by  eqns  (2.3.15)  and  (2.3.16). 


Consider  first  E and  initially  limit  attention  to  the 
m 

range  Tf» = o to  N/2  i.e.,  x = o to  L/2.  Since  E is  an  even  function  it 

in 

can  be  expanded  in  a cosine  series  of  the  form 
N/2, 


11 

S = E A cos(2mmf/N) 
m n 

n=o 


(2.3.18) 


where 


N/2 

A = (4/N)  E E cos  (2mnir/N) 
n m 


(2.3.19) 
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These  results  follow  from  eqns  (2.3.7)  and  (2.3.8)  on  replacing  N 
by  N/2.  On  substituting  for  from  (2.3.15)  in  (2.3.19)  there  results 


N/2  1 

A = (2/N)  E g cos(2mn*/N)  + 
n m 

m=o 


N/2 


+ (2/N)  I 
m=o 


11 


g„  cos(2mnir/N) 
N-m 


n = 0, • • • ,N/2  (2.3.20) 


In  the  second  series,  write 


N-m 


(2.3.21) 


The  animation  terminals  for  k are  N and  N/2,  and  since 


cos  { 2(N-k)  nir/N  } = cos(2knir/N) 


(2.3.22) 


it  can  be  seen  that  (2.3.20)  simplifies  to 


N 


11 


A = (2/N)  E g cos  (2mnir/N) 
n __  ® 

BPO 


n = 0| • • • fN/2  (2.3.23) 


Recalling  (2. 3.-11)  this  finally  becomes 


N-l 


A = (2/N)  Z g cos(2mnir/N) 


(2.3.24) 


n m 

m=o 


A similar  sort  of  calculation  shows  that,  again  for  the  half  range  m - o to 
N/2, 

N/2-1 


Z B sin  (2mnir/N) 
m , n 

n=l 


(2.3.25) 


with 


N-l 


B = (2/N)  Z g sin  (2mnu/N) 
n m=l  m 


n = l,***, N/2-1  (2.3.26) 

Thus,  using  (2.3.12),  it  follows  that  for  m = 0 to  N/2, 


N/2-1 

q = A /2  + Z {A  cos(2mmr/N)  + B sin  (2mnir/N) }+ 
m o , n n 

n=l 

(An/2/2)  cos  mir  (2.3.27) 


where  A and  B are  given  by  (2.3.24)  and  (2.3.26). 
n n 


Now  consider  the  case  of  m greater  than  N/2 
(i.e.  x in  the  range  L/2  to  L) . As  already  remarked 


E 

m 


E 


N-m 


(2.3.28) 


where  now  (N-m)  will  be  in  the  range  (0,N/2) . Thus,  one  can  use  eqn  (2.3.18) 
and  simply  replace  "m"  by  "N-m"  ? but  by  virtue  of  a relation  such  as 
(2.3.22)  it  follows  that  (2.3.18)  is  also  valid  for  the  present  range  of  m. 
Again, 


0 


m 


-0 


N-m 


(2.3.29) 


and  0 is  obtained  by  replacing  "m"  by  "N-m"  in  eqn  (2.3.25).  Since 
N-m 

sin  { 2 (N-m)  rnr/N  } = - sin  (2mnir/N)  (2.3.30) 

it  follows  that'  (2.3.25)  is  also  valid  for  m > N/2.  Hence,  the  result 
(2.3.27)  is  valid  for  all  m in  (0,N-1) . 

It  is  most  import  to  note,  though,  that  whilst  the  suffix 
m , which  is  associated  with  the  original  variable  ranges  from  0 to  N-l, 
the  suffix  n associated  with  the  transform  only  ranges  from  O to  N/2. 

The  validity  of  the  expansion  (2.3.27),  with  the  associated 
relations  (2.3.24)  and  (2.3.26)  can  be  established  directly  using 
orthogonality  relations  for  the  combined  set  of  sines  and  cosines. 


2.3.4  The  Complex  Form  of  the  DFT 

The  commonest  form  of  DFT  encountered,  and  that  which  will 
be  used  almost  universally  in  all  the  following,  is  that  obtained  from  the 
DFTT  by  replacing  the  sines  and  cosines  by  complex  exponentials?  in  fact, 
this  will  be  simply  denoted  by  "DFT". 
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The  DFTT  was  defined  by 


N/2-1 


= A /2  + E {A  cos  (2mmr/N)  + B sin(2mnir/N) } + 

o , n n 

n=l 


<V2/2) 


cos  mir 


(2.3.31) 


with 


N-l 


A = (2/N)  E g cos (2mmr/N) 
n m=0  m 


(2.3.32) 


N-l 


B = (2/N)  E g sin  (2mnir/N) 


(2.3.33) 


m=l 


On  introducing  the  complex  exponentials,  (2.3.31)  becomes 


N/2-1 


= A /2  + E C exp  (2irimn/N)  + exp  ^ + 


°'  n=l  " 


N/2-1 

+ E C*  exp  (-2nimn/N) 

n=l  n 


(2.3.34) 


where 


C = (A  - iB  ) /2 
n n n 


(2.3.35) 


C*  = (A  + iB  ) /2 
n n n 


(2.3.36) 


Consider  the  second  series.  Clearly 


exp  ( 2irim  (N-n)/N  } = exp(-2iri  mn/N) 


(2.3.37) 


and  introducing  a new  summation  index  k with 


k = N-n 

the  summation  terminals  for  k will  be  (N/2  + 1,  N-l) 


(2.3.38) 
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It  follows  from  (2.3.32)  and  (2.3.33)  that 
^N-k  = *k 


N-l 

Z C exp(2Tri  mk/N) 

k=N/2+l  * 


Since  k is  only  a dummy  suffix,  and  since  one  can  take 


B 

o 


0 


(2.3.40) 


(2.3.41) 


(2.3.42) 


(2.3.43) 


it  follows  that  eqn  (2.3.34)  can  be  re-written  as 
N-l 

q = Z C exp  (2iri  mn/N)  (2.3.44) 

m A n 
n=0 

Also,  since  eqn  (2.3.35)  now  can  be  taken  to  hold  for  all  n,  on 
substituting  from  eqns  (2.3.32)  and  (2.3.33)  into  it,  one  obtains 

N-l 

C = (1/N)  Z g exp  (-2iri  mn/N)  (2.3.45) 

n _ m 

m=0 

The  two  relations  (2.3.44)  and  (2.3.45)  constitute  a DFT  pair.  The 
second  is  sometimes  called  "the  transform"  and  the  first  "the  inversion 
formula"  in  analogy  with  the  integral  case. 
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The  validity  of  the  above  relations  can  be  established 
directly  and,  because  of  their  importance,  this  will  be  done.  Assuming 
the  expansion  (2.3.44)  exists,  multiply  both  sides  by  exp  (-2IIi  mk/N)  , 
where  k is  in  (0,N-1) , and  sum  over  m.  The  result  is 

N-l  N-l  N-l 

l g exp  (— 2rri  mk/N)  = I C { Z exp(2iri  m(n-k)/N)} 
m=0  m n=0  n m=0  (2.3.46) 

The  inner  sum  on  the  right  hand  side  of  (2.3.46)  is 
just  a geometric  series,  first  term  unity  and  common  ratio  exp  (2iri  (n-k)/N) . 


N-l 

E exp  (2irim(n-k)/N)  = {l-exp2iti  (n-k)  } / { 1-exp  (2iri  (n-k)/N) 
m=0 

(2.3.47) 


i.e. , 

N-l 

E exp  (2irim(n-k)/N)  = N 6^ 
m=0 

where  6 is  the  Kronecker  delta  defined  by 
nj> 

6=0  n / k 

nk 

= 1 n = k 

The  relation  (2.3.48)  may  be  deemed  the  analogue  of  the  Dirichlet  kernel 
for  a finite  series.  Using  (2.3.48)  in  (2.3.46)  establishes  the  required 
result  (2.3.45) . 


(2.3.48) 


(2.3.49) 


It  should  be  noted  that  this  last  proof  of  the  DFT  is 
valid  for  complex,  as  well  as  real,  gm-  Parseval's  theorem  obtained  by 
multiplying  each  side  of  (2.3.44)  by  its  conjugate  and  summing  over  m 
takes  the  form 

N-l  N-l 

lag*  = N E C C*  (2.3.50) 

ym  ym  „ n n 


Before  passing  on  attention  is  dram  to  the  fact  that 

when  g is  real  (as  was  assumed  in  the  first  derivation  of  the  DFT 
m 

given  above)  then. 


C = CJ 
n N-n 


(2.3.51) 


This  relation  is  simply  (2.3.41)  rewritten. 

2.4  A Real  DFT  as  an  Interpolation  Formula 

One  interpretation  of  a real  DFT  that  is  often  useful  is  as  a 
trigonometric  interpolation  formula.  That  this  is  a valid  interpretation 
can  be  readily  seen  as  follows  for  the  case  of  the  DFST.  (The  argument 
is  quite  analogous  for  the  other  real  DFT's). 

Consider  a function  g (x)  which  is  defined  over  the  interval 
(0,L)  by  the  relation 


where 


g(x)  = ZB  sin(nifx/L) 
n=l  n 


B = (2/N)  Z g sin  (m  mr/N) 
n m=l  m 


The  notation  g is  here  used  to  denote 
m 


g (x  ) where 
m 


(2.4.1) 


(2.4.2) 


x = m L/N 
m 


m=l, • • • , N-l 


(2.4.3) 


Now,  the  function  defined  by  (2.4.1)  exists  at  ALL  values  of  x in  (0,L) . 
Also,  it  follows  from  the  formulae  for  a DFST  that,  at  a point  x^  given 
by  eqn  (2.4.3),  the  series  (2.4.1)  will  return  the  value  g^.  Hence  the 
two  relations  (2.4.1)  and  (2.4.2)  do  indeed  provide  a formula  which 
interpolates  the  points  (*m»9m) • The  analogy  with  polynomial  interpolation 
can  be  displayed  by  substituting  from  (2.4.2)  into  (2.4.1),  and  reversing 
the  order  of  summation  when,  after  some  trigonometric  manipulations,  the 
result  becomes 


g(x)  = 


£ g. 

mr=l 


m 


S (x) 
m 


with 


(2.4.4) 


S (x)  = (1/2N)  (-l)“sin(mir/N)sin(NirxA) 

xn  

sin{n(x/L  -m/N)/2}sin{it  (x/L  + m/N)/2} 


(2.4.5) 


It  can  be  readily  established  that 


5m  <V 


S (x  ) = 1 

m m 


k / m 


(2.4.6) 


as  is  required  for  an  interpolation  function. 

Analogous  results  cam  be  established  for  the  DFCT  and  the 
DFTT;  for  example,  in  the  case  of  the  latter  it  cam  be  proved  that 
the  appropriate  interpolation  formula  is 

N-l 

g (x)  = £ g T (x)  (2.4.7) 

m m 


where 


T (x)  = (1/K)  (-l)msin(Nirx/L)  cos  {ir(x/L  -m/N)} 
m — 

sin  {x (x/L  -m/N)} 


(2.4.8) 


The  T , of  course,  have  the  properties  (2.4.6).  For  the  simple  case  of 
m 

N * 4,  they  are  shown  in  Fig.  2.4. 

It  might  be  worth  noting  that  the  above  interpolation  formulae 
can  be  used,  in  turn,  for  developing  schemes  for  numerical  integration 


and  the  like. 


- 51  - 


t 


2.5  A DFT  as  a Sampling  Formula;  Aliasing 

2.5.1  General 


In  practice  any  DFT  is  generally  used  as  an  approxim- 
ation to  some  function  defined  for  continuous  values  of  x ; if  this 
function  is  periodic  then  the  DFT  may  be  regarded  as  approximating  an 
infinite  Fourier  series  whilst  if  it  is  aperiodic  the  DFT  may  be 
regarded  as  approximating  a Fourier  integral.  Now,  a DFT  exactly 
reproduces  the  values  of  the  function  it  is  approximating  at  a finite 
number  of  equally  spaced  points;  hence,  it  can  be  regarded  as  a 
"sampling"  of  the  function  over  these  points.  In  this  Section  some 
account  is  given  of  the  errors  that  can  arise  in  the  approximation  of 
a continuously  defined  function  by  a DFT  due  to  the  necessarily  finite 
distance  between  samples. 

2.5.2  Sampling  from  a Periodic  Function 

Suppose  g(x)  is  a real,  periodic  function  defined  over 
the  basic  interval  (0,L) . Then  it  can  be  represented  exactly  by  a Fourier 
series  as  follows  : 


g(x)  = E c exp  (2irikx/L) 

k=-co  k 


(2.5.1) 


with 


c^  = (1/L)  /qL  g(x)  exp  (-2irikx/L)  dx 

and 


(2.5.2) 


c = c*  (2.5.3) 

-k  k 

the  asterisk  denoting  a complex  conjugate. 

On  the  basis  of  a subdivision  at  the  N points 


m 


mL/N 


m=0 , • • • N— 1 


(2.5.4) 


52  - 


I 

♦ 

I 


' 1 


the  DFT  representation  of  the  same  function  is 
N-l 

q = I C exp  (2irimn/N) 
m . n 
n=0 


(2.5.5) 


with 


N-l 

C = (1/N)  I g exp  (-2irimn/N) 
n _ _ m 

m=0 


(2.5.6) 


and 


C - C* 
N-n  n 


(2.5.7) 


A relation  between  the  C's  and  the  c's  can  be  obtained  as  follows.  From 
(2.5.1)  it  can  be  seen  that 


g = Z c exp  (2irikm/N) 
m . _ k 

k=-o° 


(2.5.8) 


Substituting  from  this  last  into  (2.5.6)  and  reversing  the  order  of 
summation  gives 

N-l 

C = I c {(1/N)  E exp  (2ffim(k-n)/N) } (2.5.9) 

n k=-«.k  m=0 

It  is  convenient  to  consider  the  value  of  the  inner  bracketed  sum  over 
different  ranges  of  k separately. 

k in  (0,N-1) 

From  the  basic  orthogonality  relation  (2.3.48),  the  sum  is  just 

k in  (-N,-l) 


(i) 


(ii) 
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exp  (2irim(j-N)/N)  = exp  (2irimj/N) 


the  sum  is  6^n  or,  equivalently,  6^  n_N 


(2.5.11) 


(2.5.12) 


Proceeding  in  the  same  way  for  all  such  ranges  and  inserting  the  results 
in  (2.5.9)  gives 


Cn  = Cn  + C- (N-n)  + CN+n  + c-(2N-n)  + 


(2.5.13) 


C = c + c*  + c + c*  + 
n n N-n  N+n  2N-n 


(2.5.14) 


In  order  to  interpret  this  relation  it  is  a little  simpler  to  revert  to 
real  variables.  Consider  first  the  case  when  n <_  N/2.  Then 


C = (A  - iB  )/2  ; c = (a  - ib)/2 

n n n n n n 


(2.5.15) 


where  A's  and  a's  are  the  coefficients  of  the  cosine  terms  in  the  DFTT  and 
the  Fourier  trigonometric  series  respectively;  and  likewise  the  B's  and  b's 
are  coefficients  of  the  sine  terms.  On  substituting  from  (2.5.15)  into 
(2.5.14)  and  equating  real  and  imaginary  parts  the  result  is 


An  = an  + Vn  + Vn  + a2N-n  + 


(2.5.16) 


Bn  " bn  “ bN-n  + bN+n  " b2N-n  + 


(2.5.17) 


Adopting,  for  the  moment,  the  interpolation  viewpoint  for  the  DFTT,  it  can 
be  regarded  as  representing  an  approximating  function  gg  (x)  given  by 


N/2n 

q (x)  = Z {A  cos  (2irnx/L)  + B sin  (2nnx/L) } 
s - n n 

n=0 


(2.5.18) 


On  the  other  hand,  the  exact  function  is  given  by 


g (x)  * I {a'  cos  (2irnx/L)  + b sin  (2irnxA)> 
_ n n 

n=0 


(2.5.19) 
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(In  the  above,  the  single  accent  on  the  summation  sign  indicates  that 
a weight  factor  of  1/2  is  to  be  assigned  to  the  first  term;  the  double 
accent  has  the  same  meaning  as  earlier.) 

Now  substitute  from  (2.5.16)  and  (2.5.17)  into  (2.5.18) 

and  compare  the  result  with  (2.5.19).  It  can  be  seen  that,  in  the  sample 

function,  the  coefficients  of  these  trigonometric  terms  having  a given 

value  of  n N/2)  consist  not  only  of  the  "true"  coefficients  a^  and 

b but  also  an  infinite  set  of  other  coefficients  hav  ing  suffices  out 
n 

beyond  N/2.  This  phenomenon  is  known  as  "aliasing”  or  "spectrum  folding". 
The  reason  for  the  second  name  can  be  seen  from  Fig. 2. 5.  If,  on  a strip 
of  paper,  the  number  n and  all  its  "aliases"  cure  narked  as  in  Fig. 2. 5 (a), 
and  if  then  the  strip  is  folded  at  the  points  N/2,  N,  3N/2  etc.  the 
original  number  and  all  its  aliases  sit  one  upon  the  other  as  in  Fig. 2. 5(b). 


Consider  now  the  circumstances  under  which  aliasing  does 
not  occur.  This  will  only  be  the  case  if  the  a^  and  b^  are  all  zero 
(or  ignorably  small)  for  k > N/2.  This,  then,  imposes  a constraint  on  the 
choice  of  N;  equivalently,  since  the  distance,  h,  between  adjacent  sample 
points  is  simply 


h = L/N  (2.5.20) 

it  is  a constraint  on  h.  It  is  necessary  to  choose  N sufficiently  large, 
or  h sufficiently  small,  so  that  the  aliasing  error  is  small. 

The  above  condition  is  usually  stated  in  the  context  of 
frequency  analysis.  The  last  trigonometric  term  in  (2.5.18)  is 
cos (2itNx/2L)  and  this  represents  an  oscillation  of  frequency 

fOT  = N/2L  (2.5.21) 

This  frequency  is  known  as  the  Nyquist  frequency.  The  requirement  that 
there  should  be  no  terms  in  the  series  (2.5.19)  with  n > N/2  can  be  stated 
as  follows.  "The  maximum  frequency,  F say,  in  the  true  function  must  not 
exceed  the  Nyquist  frequency  of  the  sample  function."  This  is  the  Nyquist 
criterion  for  the  prevention  of  aliasing  and  can  be  written  as 


h $ 1/ ( 2F ) 


(2.5.22) 
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In  this  form,  the  criterion  stipulates  how  closely  together  samples  must 
be  taken. 


A simple  example  of  how  aliasing  arises  is  shown  in  Fig. 2. 6. 
There  N = 4 (and  L = 1 for  simplicity) . The  function  cos  2irx  corresponds 
to  n = 1 in  the  above;  since  N-n  is  3,  the  first  alias  function  is 
cos  6irx.  As  seen  from  the  Figure,  these  two  functions  coincide  at  the 
sample  points  and  so  are  indistinguishable. 


A physical  example  of  aliasing  is  the  stroboscopic  effect. 


Finally,  it  will  be  recalled  that  in  all  the  discussion 
following  eqn  (2.5.14)  it  was  assumed  that  n ^ N/2.  However  when  n > N/2 
(but  f N-l)  it  is  easy  to  show,  by  writing  n=N-k  that  instead  of  eqn  (2.5.14) 
one  obtains  its  conjugate  complex.  The  subsequent  relations  for  the  real 
variables  (where  n is  always  less  than  N/2)  are  unaltered. 


2.5.3  Sampling  from  an  Aperiodic  Function 


It  has  just  been  shown  that  a "band-limited"  periodic  function 
(i.e.  one  whose  maximum  frequency  is  bounded)  can  be  reconstructed  exactly 
from  a number  of  equally  spaced  samples  provided  these  samples  are  sufficiently 
closely  spaced.  A similar  result  applies  to  a band  limited  aperiodic  function 
(Fig. 2. 7),  where  it  generally  goes  under  the  title  of  Nyquist's  theorem.  An 
outline  of  the  proof  is  given  below. 


If  g (x)  is  the  function,  then  the  statement  that  it  is  band- 
limited  means  that 


G(f)  = 0 


f > F 


(2.5.23) 


where  F is  considered  known.  Hence  one  can  write 


F 


g (x)  = / G (f ) exp  (2IIifx)  df 

-F 


(2.5.24) 


Suppose  now  that  G is  expanded  as  a complex  Fourier  series  over  the  basic 
interval  (-F,F) . The  result  may  be  written 


G(f ) = £ c exp  (-2irinf/2F) 

n 

n*-* 


(2.5.25) 
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with 

F 

c * (1/2F)  / G(f)  exp  (2ninf/2F)  df  (2.5.26) 

n -F 

(Here  a slightly  different  form  of  the  series  has  been  used  to  that  given 
in  Chapter  1.)  Comparing  eqns  (2.5.26)  and  2.5.24)  it  can  be  seen  that 

c = (l/2F)g (n/2F)  (2.5.28) 

n 

Consider  now  a sampling  of  g(x)  done  at  the  points  (Fig. 2. 8) 


x = n h n = 0,  1,  -1,  ••• 

where 

h = 1/2F 

and,  as  usual,  write 

g = g (nh) 
tv 

It  is  readily  established  that  the  sample  function 


(2.5.29) 


(2.5.30) 


(2.5.31) 


q (x)  = E g sine  2irF(x-nh)  (2.5.32) 

ys  n=-°°  n 

returns  the  values  gn  at  xr  provided  h is  related  to  F by  eqn  (2.5.30). 
(This  is  a valid  interpolation  formula  whatever  be  the  significance  of  F.) 

The  transform  Gg  of  gg  can  be  obtained  by  first  recalling  that  a rectangular 
pulse  of  height  1/2F,  and  extending  from  -F  to  F,  has  the  transform 
sine  2irFx;  and,  so,  vice-versa.  In  order  to  obtain  the  transform  of  the 
translated  sine  functions  occurring  in  eqn  (2.5.32)  the  shift  theorem  is  then 
used.  In  this  way  it  is  established  that 


G (f)  = (1/2F)  ? g exp  (-2irifnh) , |f|  < F 

s n=-°°  n 

=0  |f|  > F 


(2.5.33) 


As  long  as  h is  chosen  in  accord  with  eqn  (2.5.30),  it  can  be  seen  that 
G coincides  with  G.  Since  a function  is  uniquely  determined  by  its 


57 


transform  (through  the  inversion  integral)  it  follows  that  the  sample 
function  coincides  with  the  original  one. 

If  h is  chosen  smaller  than  is  required  by  (2.5.30) 
the  argument  remains  valid.  A smaller  value  of  h may  be  considered 
to  correspond  to  a larger  value  F*  of  the  bandwidth;  however,  a function 
which  has  a bandwidth  F can  certainly  be  deemed  to  have  a bandwidth 
F1  > F. 

On  the  other  hand  if  h is  chosen  larger  than  required 
by  (2.5.30)  the  argument  breaks  down.  The  transforms  G and  G will 
no  longer  coincide,  and  it  can  be  shown  that  the  errors  involved  in  using 
Gg  can  be  interpreted  in  terms  of  aliasing. 

2.6  Truncation  and  End  Effects;  Windows 

2.6.1  General 

Consider  an  aperiodic  function  g(x)  which,  by  definition, 
is  given  over  an  infinite  range  of  x.  When  approximating  this  function  by 
a DFT,  one  necessarily  is  working  only  with  some  finite  range  (0,L)  of  x. 
The  question  is,  what  error  is  involved  in  this  "truncation"  of  the  true 
function.  This  will  now  be  taken  up. 

It  turns  out  that  the  same  sort  of  error  arises  when  the 
DFT  is  used  to  approximate  a periodic  function  but  when  the  length  L does 
not  coincide  with  the  period  (or  seme  integral  multiple  of  the  period.) 
This,  of  course,  might  have  been  expected  on  the  basis  of  the  remark  made 
earlier  that  an  aperiodic  function  can  be  thought  of  as  a periodic  one 
but  with  infinite  period.  ' 

2.6.2  Truncation  of  Interval 

If  the  true  function  under  consideration  is  g(x)  then 
the  truncated  function  gfc (x)  is  given  by  (Fig.  2.9) 

gt(x)  = g(x)  o < x < L 


0 


elsewhere 


(2.6.1) 


I 


It  is  clear  that  one  can  write 

9t(x)  * 900  w(x)  (2.6.2) 

where  w(x)  is  the  "rectangular  window"  function  defined  by 


w (x)  = 1 


o < x < L 


= 0 


elsewhere 


(2.6.3) 


(The  terminology,  of  course,  arises  from  the  fact  that  the  truncated 
function  can  be  considered  as  a restricted  "view"  of  the  true  one  through 
the  "window" . ) 

Since  (2.6.2)  is  a product,  the  transform  Gfc  can  be 
written  as  a convolution.  The  transform  of  w was  shown  in  Section  1.5.2 
to  be 


W(f)  = L exp  (-irifL)  sinc(irfL)  (2.6.4) 

Hence,  by  the  convolution  theorem, 

v 

G (f)  = /”  G(u)  L exp  l-iri  (f-u)L)  sine  $1T(f-u)L}  du  (2.6.5) 

t —00  * N 

In  order  to  interpret  this  result,  consider  the  case  where 
g(x)  is  a pure  (complex)  oscillation  corresponding  to  a frequency  f i , say 
x . e. , 

g(x)  = exp  (2trif1x)  (2.6.6) 

Then, 

G(f)  = 6 (f-f i)  (2.6.7) 

and 

Gt(f)  = L exp  (-iri(f-fi)}  sine  {ir(f-fi)L)  (2.6.8) 

In  particular,  the  magnitude  of  Gfc  is  given  by 

| G (f ) | = L | sine  {n(f-fi)  L}  | (2.6.9) 


T-^ 
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Comparing  (2.6.7)  and  (2.6.9)  it  can  be  seen  that,  whereas  the  transform 
of  the  true  function  is  a "spike"  at  the  point  i\,  the  transform  of  the 
truncated  function  gives  a "diffraction  pattern"  centred  about  fj 
(Fig. 2. 10).  If  one  were  dealing  with  a function  that  comprised  the  sum 
of  two  oscillations  whose  frequencies  were  close  together  then  the  two 
diffraction  patterns  would  overlap;  in  such  a case  it  may  not  be  clear 
from  a study  of  the  transform  what  the  true  situation  is.  This 
phenomenon,  where  truncation  effects  lead  to  the  replacement  of  a "spike" 
by  a diffraction  pattern,  is  known  as  "leakage". 

2.6.3  DFT  of  a Function  "Not  Periodic  in  the  Window" 


A DFT  defined  over  the  interval  (0,L)  represents,  outside 
that  interval,  the  periodic  repetitions  of  the  original  function.  Here 
one  is  concerned  with  the  case  where  the  original  function  is  periodic  but 
its  period  (or  an  integral  multiple  of  it)  may  not  be  equal  to  L.  Simple 
examples  of  this  situation  are  shown  in  Fig. 2. 11;  there,  the  first  example 
sin  (2irx/L)  does  have  a period  L but  the  second  sin  (1.5irx/L)  does  not. 
In  the  first  case  the  extended  function  is  continuous  whilst  in  the  second 
the  extended  function  has  discontinuities  at  the  end  points.  It  is  these 
discontinuities  which  cause  trouble. 


-Consider  the  complex  oscillation  of  frequency  fj  (and  hence 
period  Lj  = 1/fi)  given  by 


g(x)  = exp  (2irif1x) 


(2.6.10) 


Recalling  the  general  formulae  for  a DFT,  namely, 


g = I C exp  (2irimn/N) 
m _ n 

n=0 


(2.6.11) 


N-l 

C = (1/N)  I 

^ r\ 


N-l 

I g exp  (-2irimn/N) 
_ m 
m=0 


(2.6.12) 


it  is  easily  seen  that  if 


■ k/L  (k  integer) 


(2.6.13) 


then 

C =1  , C = 0 n * k 

K n 


(2.6.14) 


However , suppose 

fx  * (k  + p)/L  , o < p < 1 (2.6.15) 

Now 


g = exp  {2iri(k+p)m/N  } (2.6.16) 

m 

On  substituting  from  (2.6.16)  into  (2.6.12),  and  noting  that  the  result 
is  a geometric  series  similar  to  eqn  (2.3.46),  it  is  found  that 

Cn  = (1/N)exp{iri (k+p-n)  (1-1/N)  }sin{Tt^t+p-n)  }/sin{ir (k+p-n)/N> 

(2.6.17) 

Confining  attention  to  the  magnitude  of  the  coefficients,  it  can  be  seen 
that 


Ic^l  = (1/N)  | sin{ir  (k+p-n)  } /sin{ir  (k+p-n)/N)  | (2.6.18) 

The  ratio  of  the  two  sines  occurring  in  (2.6.18)  has  a behaviour  generally 
similar  to  the  sine  function.  In  general  all  the  are  non-zero, 

although  the  largest  coefficients  are  and  C^+1  as  might  be  expected. 


( 


<■ 


* 


- _ . ^ H -WWIWBijmiiijmjt.n.i 

- 61  - 

Again,  this  represents  "leakage"  compared  with  the  "true"  situation  as 
represented  by  something  like  (2.6.14).  Leakage  may  be  considered  here 
as  the  error  arising  from  an  incorrect  choice  of  the  length  of  the 
basic  interval  or  because  "the  function  is  not  periodic  in  the  window". 
The  evaluation  of  eqn  (2.16.18)  for  the  case  N=16,  k+p=3.5  is  shown  in 
Pig. 2.12. 


2.6.4  Reduction  of  Leakage  by  Hanning 

As  already  seen,  the  problem  of  "leakage"  is  associated 
with  end-effects  which  arise  when  the  function  is  not  periodic  in  the 
window.  One  way  of  ameliorating  this  problem  is  to  use  the  concept  of 
hanning  already  introduced  in  the  Appendix  to  Chapter  1.  It  will  be 
recalled  that  there,  instead  of  dealing  with  the  original  function  g(x), 
one  considered  a modified  function  p(x)  say,  given  by 

p(x)  = g(x)  (1-cos2itx/L)/2  (2.6.19) 

The  modified  function  will  be  periodic  in  the  window.  The  discrete  form 
of  eqn  (2.6.19)  is,  of  course, 

p = g (1-cos  2irm/N)/2  (2.6.20) 

m m 

By  writing  the  cosine  in  terms  of  exponentials  and  applying  the  general 

fo  -mula  for  a DFT  to  p it  is  raadily  found  that,  if  C denotes  the 

m n 

transform  of  g,  and  C1^  denotes  the  transform  of  p,  then 

C1  = -C  /4  + C /2  - C ,/4  (2.6.21) 

n n-1  n n+1 

where  n is  in  the  range  (0,N-1) . (Values  of  C for  subscripts  lying 
outside  this  range  are  obtained  from  periodicity  considerations) . The 
results  of  applying  eqn  (2.6.21)  to  eqn  (2.6.17)  are  shown  in  Fig. 2. 12 
and  it  can  be  seen  that  the  leakage  has  been  decreased.  As  with  the 
continuous  case,  if  the  modified  function  is  reconstructed  using  the 
hanned  coefficients,  ♦■-he  original  function  can  be  obtained  from  eqn  (2.6.20) 


Some  light  is  thrown  on  hanning  by  considering  the 


continuous  transform  of  the  "Hann  window" , 


w(x)  a (1-COS  21f  x/L)/2 


o ^ x i L 


X > L 


(2.6.22) 


After  some  routine  calculations  it  is  found  that 

W(f)  = { (-1)  (L/2)exp(-irifl)sin(irfL)  } / {irL3f  (f-l/L)  (f+l/L)  } 

(2.6.23) 

The  modulus  of  this  function  has  been  shown  in  Fig. 2. 13  along  with  the 
corresponding  result  for  a rectangular  window  (where  the  transform,  of 
course,  was  basically  sincivfL) . The  main  points  are 

(i)  the  height  of  the  main  lobe  is  here  half  that  for 

the  rectangular  window. 


the  width  of  the  main  lobe  is  here  twice  that  for 
the  rectangular  window  (which  is  a detrimental 
effect  as  regards  the  resolution  of  close  frequencies) , 


the  maximum  amplitudes  in  the  minor  lobes  drop  off 
much  more  rapidly  than  for  a rectangular  window 
and  this  reduces  leakage.  (Basically,  these 
amplitudes  are  now  dropping  of  as  1/f3  rather  than  1/f 
for  the  rectangular  window.) 


2.7.  The  Discrete  Convolution 


2.7.1 


The  Convolution  Theorem 


It  is  possible  to  derive  a convolution  theorem  for  a DFT 

which  is  quite  analogous  to  that  for  a continuous  transform.  Suppose  that 

it  is  required  to  evaluate  the  transform  B^  of  a point  function 

which  can  be  considered  as  the  product  of  two  other  point  functions,  h 

in 

and  y , say,  whose  transforms  H , Y are  known;  the  question  is,  to 
m n n 

express  B in  terms  of  H and  Y . 

n n n 
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It  is  given  that 


z = y h 
m ^m  m 


m = 0, • • • , N-l 


(2.7.1) 


where 


m 


m 


E Y^  exp (2-rrimj/N) 


(2.7.2) 


£ 


exp(2nimk/N) 


(For  brevity,  the  summation  terminals  of  0 and  N-l  have  been  omitted.)  From 
the  definition,  and  using  (2.7.1), 


8 = (1/N)  E y h exp  (-2irimn/N) 

n m m m 


(2.7.3) 


On  substituting  from  (2.7.2)  into  (2.7.3),  and  altering  the  order  of 
sumnation,  one  gets 


8 = (1/N)  Z E Y H {Z  exp(2irim(j+k-n)/N) } 


m 


(2.7.4) 


j k j k 

From  the  orthogonality  relation,  the  inner  sum  has  the  value 


N <S, 


Hence, 


k,n-j 


N-l 


8 = . E Y.H 

n j«o  j n-] 


(2.7.5) 


This  last  summation  is  termed  a discrete  convolution  of  the  function*  H 
and  Y. 


If,  instead  of  (2.7.1),  one  is  given  that 


8 = Y H 

n n n 


(2.7.6) 


then  a parallel  argument  shows  that 


N-l 

z = (1/N)  E y h . 
m k=o  k m-k 


(2.7.7) 
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2.7.2  Wrap-Around  Errors 


Consider  again  the  discrete  convolution  given  by 


N-l 

Zm  = klo  yk  Vk 


m=0, • • • ,N-1 


(2.7.8) 


where,  for  brevity,  the  constant  factor  (1/N)  has  been  omitted.  First 
note  that  the  suffix  (m-k)  takes  values  in  the  range  -(N-l)  to  (N-l). 
Secondly,  in  all  calculations  using  the  DFT  it  is  implicit  that  any 
function  appearing  is  of  period  N.  In  the  present  case,  then 


h . = h„  . 
-3  N-J 


(2.7.9) 


These  two  points  create  difficulties  when  a discrete  convolution  is  used 
as  an  approximation  to  the  continuous  convolution  integral. 

The  continuous  convolution  to  which  (2.7.8)  may  be  regarded 
as  some  sort  of  approximation  will  be  of  the  form 


z(x)  = / y(C)  h(x-C)  <3? 


(2.7.10) 


If  this  is  approximated  by  a sum  using  the  same  sub-division  as  employed 
in  the  discrete  case,  but  extended  at  each  end,  this  will  be  of  the  form 


z 

m 


m-k 


m=-». 


(2.7.11) 


where,  again  in  the  interests  of  brevity,  a constant  factor  giving  the 
length  of  sub- interval  has  been  omitted. 

If  the  infinite  sum  represented  by  (2.7.11)  is  to  coincide 
with  the  finite  one  (2.7.8)  for  values  of  m in  (0,N-1)  i.e.,  "in  the 
window" , then  clearly  at  least  one  of  the  functions  appearing  in  the 
convolution  must  be  zero  over  "most  of  the  range".  From  here  on,  it  is 
simplest  to  consider  a concrete  case.  Hence,  take  N=8  and  suppose  that 
whilst  y is  in  general  non-zero  everywhere,  h is  only  non-zero  at  the 
points  0,-1, -2  and  -3  (i.e.,  only  N/2  values  of  h are  non-zero). 
Considering  the  infinite  sum  first,  the  only  non-zero  terms  occurring  are 
those  denoted  by  an  asterisk  in  Fig. 2. 14 (a).  For  example 
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z = yh  + y h +yh  + yh  (2.7.8) 

2 20  3-1  4-2  5-3 

and 

z =yh+yh  +yh  +yh  (2.7.9) 

6 607-1  8~2  9 “3 

Turning  now  to  the  discrete  case,  an  analogous  table  is  shown  in  Fig. 2. 14(b). 

Comparison  of  the  two  tables  shows  that  for  zq  to  z^  inclusive,  the  infinite 

and  finite  sums  provide  the  same  result.  On  the  other  hand,  for  z^  to  z^ 

the  two  sums  give  different  results.  For  example  in  z , as  seen  from 

6 

eqn  (2.7.9),  the  infinite  sum  has  terms  involving  y and  y whilst  the 

8 9 

finite  sum  has  no  such  terms.  However,  the  finite  sum  as  well  as  returning 

the  first  two  terms  in  (2.7.9),  returns  terms  in  y and  y since  it  assumes 

° 1 

that 


h = h , h = h (2.7.10) 

6-25-3 

I 

and  these  are  non-zero.  This  type  of  error  is  known  variously  as  "wrap- 
around" or  "overlap"  or,  sometimes,  as  aun  "end-effect".  Because  of  it, 
valid  convolution  results  are  only  obtained  for  "half  of  the  window".  If 
one  wished  to  obtain  the  missing  (or,  rather,  presently  invalid)  results 
it  would  be  necessary  to  repeat  the  calculation  with  the  window  differently 
located . 

It  might  be  noted  that  if  a different  number  of  non-zero 
values  of  h had  been  used,  a different  number  of  valid  convolution 
results  would  have  been  obtained.  (Increasing  the  number  of  non-zero  values 
of  h decreases  the  number  of  valid  results;  in  practice  it  is  usual  to 
take  N/2  values  as  above) . A detailed  account  of  convolution  calculations 
has  been  given  by  Brigham  (6) . 
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FIG.  2.3  EVEN  AND  ODD  POINT  FUNCTIONS  FOR 

I g(x)  * exp  (x/L)  - 1,  0<x<L  WITH  N = 10 

g(x  + L)  - g(x) 


(a)  Aliases  of  n 


(b)  Spectrum  folding 


FIG.  2.5  ALIASING  OR  SPECTRUM  FOLDING 


0 sample  point 


FIG.  2.6  ALIASING  OF  TWO  COSINE  FUNCTIONS 


FIG.  2.13  TRANSFORMS  OF  RECTANGULAR  AND  HANN  WINDOWS 


THE  FAST  FOURIER  TRANSFORM  ALGORITHM 


N-l 

G = G ( f ) = (1/N)  E g exp  (-2irimn/N) 
n n m 

m=o 


with  an  inverse  given  by 


N-l 

g = g(x  ) = E G exp  (+2irimn/N) 
mm  n 

n=o 


Using  the  common  notation 


exp  (-2Tri/N)  = cos(27r/N)-i  sin(2ir/N) 


for  the  N th  root  of  unity  eqns  (3.1.1)  and  (3.1.2)  may  be  written 
concisely  as  follows  : 


Clearly  any  method  for  calculating  the  transform  may  easily  be 
adapted  to  calculate  the  inverse  by  using  the  complex  conjugate  of  the 
exponential  terms  W™  and  deleting  the  factor  1/N.  In  the  following 
discussion  of  computational  procedures  it  is  therefore  sufficient  to 
consider  only  the  forward  transform  without  the  factor  1/N  i.e., 


In  many  cases  it  is  convenient  to  write  the  N equations  (3.1.6)  in 
matrix  notation  : 


g = Mn  g 

where 


“»■ 


"w° 

w° 

w° 

...  w° 

w° 

w1 

w2 

...  w”"1 

w° 

w2 

w4 

...  w2(N 

. 

m 

• • • , 

w°  «N-1w2<N-1)... 


(3.1.7) 

iT 


(3.1.8) 


The  matrix  representation  shows  that  the  direct  calculation  of  the  DFT 
(or  its  inverse)  requires  the  multiplication  of  the  N-element  complex 
vector  g by  the  N x N complex  matrix  to  form  the  N complex 
elements  of  the  transform  vector  G.  Assuming  that  the  N2  elements  of 
the  matrix  Mn  have  been  calculated  once  and  for  all,  the  complete 
matrix  multiplication  requires  N2  operations  of  complex  multiplication 
and  addition.  In  typical  practical  problems,  N may  be  of  the  order  of 
103  and  so  N2  is  of  the  order  of  106.  In  the  past,  the  usefulness  of 
the  DFT  as  a numerical  tool  was  limited  by  the  fact  that  the  confutation 
time  was  proportional  to  N2. 

The  publication  of  the  paper  by  Cooley  and  Tukey  (1)  revolutionised 

the  approach  to  wave-form  analysis  (i.e.  DFT  calculations) . They  showed 

that  if  the  number  of  samples,  N,  is  a conposite  number  having  factors 

r r , •••  r.  , the  number  of  operations  required  to  calculate  the  DFT 
1 2 * 
may  be  reduced  from 
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Their  method,  and  similar  later  ones,  for  calculating  the  DFT  have 
become  known  as  the  Fast  Fourier  Transform  (FFT) . It  should  be 
emphasized  again  that  the  FFT  is  simply  an  algorithm  for  calculating 
the  DFT;  the  numerical  result  is  the  same  as  if  the  full  matrix  multi- 
plication had  been  performed. 

The  advantage  of  the  FFT  is  a dramatic  reduction  in  computational 
time  achieved  by  use  of  an  algorithm  which  is  an  efficient  rearrangement 
of  the  process  used  for  matrix  multiplication.  A second,  but  less 
important  advantage,  is  a possible  saving  in  computer  memory  by  the  use 
of  "in-place"  computation.  This  last  aspect  will  be  discussed  in  a 
subsequent  chapter. 

Here,  the  basic  two  factor  algorithm  will  be  derived  using 
summation  notation.  Most  references  use  this  notation,  presumably  because 
of  the  origin  of  the  DFT  as  an  approximation  to  the  classical  Fourier 
integral.  However,  as  the  matrix  notation  has  special  advantages,  the  FFT 
algorithm  will  be  described  in  this  form  as  well.  Finally,  a description 
of  multifactor  FFT's  will  be  given. 


The  Basic  Fast  Fourier  Transform  Algorithm 


The  FFT  is  an  algorithm  for  computing  the  N elements  G(n)  in 
(3.1.6)  using  less  than  N2  operations.  The  method  is  only  applicable 
when  N is  non-prime. 

The  algorithm  for  N even,  which  exploits  the  factor  2,  is  often 
given  as  an  introductory  example,  e.g.  refs. (2)  and  (3),  since  it  may 
readily  be  derived  from  first  principles.  As  this  approach  rather  obscures 
the  generality  of  the  FFT,  the  algorithm  derived  below  is  for  N having 
two  arbitrary  factors. 

If  N = r . r then  the  N elements  g 

1 2 m 

of  the  original  function  and  the  N elements  G^  of  the  transform  may 


♦ elements.  The  basis  of  the  FFT  algorithm  is  to  group  the  elements  of  the 

original  function  by  the  first  method  and  the  elements  of  the  transform 


Ml 


by  the  second.  It  will  be  convenient  to  consider  the  successive  groups 

of  the  elements  g as  being  the  columns  of  an  r x r.  array  and  the 
m * > 

successive  groups  of  the  elements  as  the  columns  of  an  t ^ x r^ 


array. 


With  this  convention  the  index  m (m  = 0,  1,  •••  N-l)  which  tags 


the  original  elements  is  written 


m = m + m r 
o i 2 


where  "row"  index  m = 0,  1,  ...  r2  - 1 

o 


and  "column"  index  m = 0,  1,  •••  r -1 

1 1 


Hence  one  can  write 


g(m)  - g(m  , m ) 
1 o 


(3.2.1) 


Analogously,  the  index  n which  tags  the  elements  of  the  transform  is 


written 


n = n + n r 
o l 1 


where  "column"  index  n = 0,  1,  — r -1 

o 1 


(3.2.2) 


and  "row"  index 


n = 0,  1,  * • • r -1 
1 2 


Thus,  G(n)  = G(n^,  nQ) 


For  example,  if 


N = 12  = 4 X 3,  the  two  groupings  are  as  shown  in  Fig.  3.1. 

Using  the  above  representations  of  the  indices  m and  n,  the 
exponent  n m in  eqn  (3.1.6)  becomes 


I 


i 
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n.m  = (n  + n,r,)  (m  + m.r.,) 
o ll  o i 2 


r,n,m  + n m + r_n  m,  + r.r_n.m. 
1 1 o oo  2 o 1 1211 


Thus,  the  factor  W„  of  eqn  (3.1.6)  becomes 
N 


nm  r,n,m  n m r_n  m, 
,,  1 1 ° „ o o „ 2 o 1 

WN  * WN  WN  WN 


(3.2.3) 


since 


rlr2nlml  . N nlml 
WN  = (WN  > = 1 


With  the  above  notation  the  expression  (3.1.6)  for  the  DFT  becomes 


r -1  r -1  r.n.m  n m r n m 
2 1 1 1 o oo  2ol 


G(n  ,n  ) 
1 o 


z „ z w„ 

m =0  m, =0  N 
o 1 


\ WN 


g(m  ,m  ) 
1 o 


(3.2.4) 


with  n = 0,  1,  ...  r2“l 


n = 0,  1,  . . . r -1 
o 1 


This  double  sum  (3.2.4)  is  the  basis  of  the  FFT.  However,  there  are  two 
similar  but  quite  distinct  forms  that  the  sum  may  take. 


Firstly,  using  the  fact  that 


exp  (-2nir2/r1r2)  = exp  (-27ri/r1)=Wr 


the  double  sum  may  be  reorganised  as 


r —1  r —1 

2 (n  +n  r )m  1 n m 

G(n.,n  )=  ° 11  ° { Z W ° 1 g(m  ,m  )} 

1 o m =0  N m,=0  r.  1 o 

o 11 


with  n^  = 0,  1,  ...  r2“l 


(3.2.5) 


. ipimiwu.uiHwi  .1.111.  Mimiuiimmu 
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nQ  * 0,  1, 


rl  -1 


This  relation  is  the  essence  of  the  Cooley-Tukey  version  of  the  FFT 
algorithm,  also  referred  to  as  "decimation  in  time".  The  inner  sums 


which  depend  only  on  mQ  and  nQ,  and  are  performed  first  are  : 


a (n  ,m  ) 
o o 


r,-l 

1 n m 

£ n W ° g (m  m ) 
m^-u  r^  1,  o 


(3.2.6) 


with 


nQ  = 0,1, •••  rx-l 


Each  of  these  expressions  is  precisely  an  r -point  DFT  of  each  of  the 

1 2 

r2  rows  of  g in  Fig. 3.1.  Since  each  rn -point  DFT  requires  r 


i 1 

operations,  the  total  number  of  operations  required  to  form  all  the  inner 


sums  is  r_.r,  = Nr,.  For  each  value  of  n = n + n,r,  in  (3.2.5)  the 
i l l oil 


outer  sum  is  : 


r2_1 


G(n,,n  ) = E W^o  a (n  ,m  ) 
1 o m =0  N o o 


(3.2.7) 


with 


n = 0,  1 ...  N-l 


Each  sum  is  a weighted  sum  of  r2  of  the  N results  so  far  formed.  The 


outer  suir^  thus  involve  r2  operations  carried  out  N times  (one  for  each 


value  of  . The  grand  total  of  operations  is  therefore  only  N(r  + r„) 

2 1 
rather  than  N = N r^r2  for  the  direct  evaluation  of  the  N-point  DFT. 


The  two  operations  specified  by  (3.2.6)  and  (3.2.7)  may  be  described 
by  the  use  of  a signal  flow  graph.  The  Cooley-Tukey  algorithm  for 
N = 12  = 4 X 3 is  shown  in  Fig. 3. 2. 


The  alternative  rearrangement  of  the  double  sum  (3.2.4)  is 


r2~1 


fit  n 1 (m  +m.r1  )n 

G(n.,n  ) = En  W ° 1 {I  WN  ° g(m.,m  )} 

1 o m =o  r in.  =0  N ± o 

O Z 1 


(3.2.8) 


/ 


with  n.  = 0,1,  ••*,  r2  -1 


This  relation  is  the  essence  of  the  Sande-Tukey  version  of  the  FFT 
algorithm,  also  referred  to  as  "decimation  in  frequency"  (4) . In  this 
version  the  weighted  sums  of  the  N elements  of  g taken  r^  at  a time 
are  formed  first  as  specified  by  the  inner  bracket.  This  requires  Nr^ 
operations.  The  outer  sums  are  then  r2 -point  DFT's  of  which  there  are 
r^  in  all.  The  grand  total  of  operations  is  again  N(r^+r2).  The  Sande- 
Tukey  algorithm  for  N=12=3X4  is  shown  in  the  signal  flow  graph  of  Fig. 3. 3. 

At  first  sight  the  two  versions  may  appear  to  differ  merely  in  the 
interchange  of  the  roles  of  r^  and  r^.  However,  also  involved  is  an 
interchange  of  the  roles  of  m and  n.  In  fact,  the  Sande-Tukey  version  is 
the  inverse  of  the  Cooley-Tukey  version.  This  inverse  relationship  is  the 
reasoning  behind  the  terms  "decimation  in  time"  and  "decimation  in 
frequency".  In  the  Cooley-Tukey  version,  the  r^-point  DFT's  of  the  inner 
sum  are  performed  on  elements  of  the  original  or  "time"  function  spaced 
r2  apart,  whereas  in  the  Sande-Tukey  version  the  results  of  the  final 
point  transforms  generate  elements  of  the  transformed  or  "frequency"  function 
spaced  r^  apart.  This  inverse  relationship  will  be  discussed  further  in  a 
later  section.  However,  at  this  stage  it  is  important  to  emphasize  that,  in 
the  present  context,  the  two  versions  are  equally  efficient. 

3.3  Matrix  Description  of  FFT  Algorithm 

In  the  previous  section  the  basic  two-factor  algorithm  was  derived 
using  summation  notation.  However  matrix  notation,  although  not  used  by 
most  authors  in  this  field,  has  several  advantages.  For  example,  any 
operations  involved  are  more  apparent  and  manipulation  of  the  matrix  factors 
is  easy.  For  this  reason  a matrix  description  of  the  FFT  algorithm  will  be 
given  below. 

The  summation  forms  of  the  FFT  developed  in  Section  3.2  may  be 
concisely  expressed  as  a factorisation  of  the  matrix  The  elements  of 
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%■ 


» 


(At  this  stage  note  that 

(i)  within  each  partition,  the  first  factor  of  each 
element  is  the  same, 

(ii)  the  last  factor  in  each  element  is  the  same  at 
corresponding  positions  in  each  partition, 

(iii)  the  middle  factor  of  each  element  is  the  same  in 
each  row  within  each  partition,  and  is  also  the 
same  as  those  in  the  corresponding  rows  of  all 
other  partitions  in  the  same  column.) 


e)  factorise  the  individual  partitions  and  hence  the  complete 

matrix.  (This  operation  is  possible  because  of  the  repetitions 
specified  above) . 

D 

The  ri^ht  hand  factor  is  a replication  of  the  matrix  Mj  along  the 

diagonal  partitions.  Normally  M.  would  be  written  in  terms  of  W but  here 

3 Z 3 

it  is  expressed  in  terms  of  W&  where  w^=w6  • Consequently  each  element  has 
been  raised  to  the  secsnA  power  i.e.,  its  exponent  has  been  multiplied  by  7.. 

The  middle  factor  is  a diagonal  matrix  whose  (3m  + n ) th  element  is 

o o 

given  by  W,nomo  where  n = 0,1,2  and  m =0,1 

DO  O 

I •] 

The  l«£t  hand  factor  contains  the  elements  of  the  matrix  M repeated 

* 

along  the  diagonal  in  each  partition.  Again  would  normally  be  written  in 
terms  of  but,  analogously  to  the  first  case,  each  element  is  raised  to  the 
tV>ird  power  since  it  is  written  in  terms  of  W . 

D 


Pease  (5)  has  suggested  the  use  of  the  Kronecker,  or  direct,  products 
of  matrices  to  represent  the  replications  shown  in  the  right  and  left  factors. 
The  Kronecker  product  of  a matrix  A and  the  3x3  matrix  B is 


'bll 

A 

b^  A 

b13 

A 

(3.3.2) 

b21 

A 

b22  A 

b23 

A 

b31 

A 

b32  A 

b33 

A 
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With  this  notation  the  two  forms  of  replication  expressed  by  the  right  and 
left  factors  in  Fig. 3. 5 may  be  simply  expressed.  The  left  factor  is 
I3  x M2  and  the  right  factor  is  M3  x I2<  Thus  the  Cooley-Tukey  version  of 
the  6-point  transform  may  be  written  as 


G * Mg  g 


d3  x m2)  d3  2(m3  x i2>  p3  2 9 


(3.3.3) 


The  general  case  where  N = rn . r_  may  be  written  with  the  same  notation  as 


G = M g 
r r 
12 


= (I  X M ) D 


(M  x I ) g 


rl  r2  rl'r2  rl  r2  rl'r2 


(3.3.4) 


where  the  (n  + m r,)th  element  of  the  diagonal  matrix  D is  given  by 

O o 1 


W^°n°  where 


n = 0, 1, • • • r.  -1 
o J. 


m = 0,1, • • • r_  -1 
o * 


The  ahnvp  matrix  description  has  been  developed  for  the  Cooley-Tukey 
FFT.  The  Sande-Tukey  version  may  be  derived  in  a similar  manner  but  in 
that  version  the  rows  rather  than  the  columns  of  are  permuted  by  pre- 
multiplication by  P_1  = P . This  has  the  effect  of  producing  the 

rir2  r2rl 

same  matrix  factors  but  written  in  the  reverse  order,  viz. 


G = P„  _ (Mx  I ) D (I  x M ) g 

2 1 1 2 rl'  2 1 2 


The  inverse  of  (3.3.5)  is  given  by 


(3.3.5) 


g = (I 


-1  „ -1 


-1  -1  ± 


x Mr  ) D * (M  XI)  P G 

1 2 1'  2 1 2 2'  1 


Using  the  property  of  the  Kronecker  products  that 


(A  x B)"1  = A-1  x B-1 


/ 
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this  last  becomes 


g = (I  x M_1)  D_1  (M-1  x I ) P G 

rl  r2  rl,r2  rl  r2  l'r2 


(3.3.6) 


When  this  equation  is  compared  with  the  Cooley-Tukey  version  (3.3.4)  the 
terms  "decimation  in  time"  and  "decimation  in  frequency"  become  clear. 

In  this  description  of  the  FFT  algorithm,  matrix  notation  has 
been  used.  As  has  been  demonstrated,  the  advantages  of  this  notation 
are  that  it  is  compact  and  that  it  is  completely  general.  Furthermore, 
it  allows  reorganisations  of  the  basic  algorithm  to  be  developed  by  making 
use  of  the  properties  of  ordinary  and  Kronecker  matrices. 

In  this  section  the  basic  Cooley-Tukey  algorithm  has  been  described 
in  matrix  notation  using  the  Kronecker  matrix  product.  This  is  equivalent 
to  the  more  commonly  quoted  summation  form.  The  matrix  approach  has  also 
been  used  to  reveal  the  inverse  relationship  between  the  Cooley-Tukey  and 
Sande-Tukey  versions  of  the  FFT  algorithm. 

Only  the  two-factor  algorithm  has  been  described  at  this  stage.  This 
is  fundamental  to  the  more  useful  multi-factor  FFT  which  is  treated  in  the 
next  Section.  The  matrix  approach  will  prove  especially  advantageous  there, 
allowing  the  iterative  development  for  more  than  two  factors  to  be  written 
explicitly. 

3.4  Multi-factor  FFT  Algorithm 

The  use  of  the  matrix  notation  facilitates  the  development  of 
algorithms  when  N has  many  factors.  In  this  Section  an  expression  is  derived 
for  the  multi-factor  Cooley-Tukey  algorithm.  This  expression  is  then  used  to 
derive  the  algorithm  when  all  the  factors  are  equal  to  a general  radix  r. 

The  two  most  important  cases  of  radix-2  and  radix-4  are  then  discussed. 


In  what  follows 
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It  is  therefore  more  convenient  to  consider  the  initial  factorisation 
in  terms  of  the  factors  r1  and  N/r^  With  this  change  in  notation  the 
two-factor  algorithm  (3.3.4)  may  be  rewritten  by  direct  substitution  as: 


G " (IN/r1  x V Vr^r/Vr,  X 

Since  N/r^  is  also  non-prime  it  may  be  written  as 

N/r^  = (r3.r4. . .r^) .r2  = (N/fr^.r^  } .r2 


(3.4.1) 


The  factor 
and  becomes 


in  (3.4.1)  may  thus  be  factorised  in  a similar  manner 


Vr1  (1N/r1r2x  “r^  ^ V^r 


The  Kronecker  product  has  the  following  properties 


(A  B C)  x I = (A  xl)  (B  x I)  (C  x I) 


2,r2 


(3.4.2) 


and 


1 x I.  = I , 
a b ab 


Using  these  with  eqn  (3.4.2)  it  follows  that 


Vx  I = (I  , xM  xl)  (D  . *■  i 

rl  rl  N/^rir2  r2  rl  1 2'  2 1 

(M  . x I ) (PM/,  _ _ x I ) (3.4.3) 

N/rir2  rlr2  N/rir2'r2  rl 

When  the  expression  (3.4.3)  is  substituted  into  (3.4.1) the  FFT  becomes 


x I ) . 


G = { I., . x M ) D . 

N/rx  r1  N/rj/ri 


(I„  / x M x I ) (D  . x I ) 

N/rlr2  r2  rl  N/rlVr2  rl 


I 


<Vrir2  X V2) ' 


^PN/r^r2»r2  x Xr^  PN/r^,r  9 


(3.4.4) 


This  process  of  substitution  may  be  continued  until  all  the  factors  of 
N (and  the  authors)  are  exhausted.  As  an  example,  take 

N=  12  =3x2x2 

The  expression  (3.4.4)  then  becomes 


G = (I4  x M3)  D43. 


(I2  x M2  x I3) (D2  2 x I3) , 


(M2  x I6) . 


^P2,2  X *3^  P4,3  g 


(3.4.5) 


If  all  the  factors  of  N are  equal  the  complete  factorisation  is 
known  as  the  radix-r  FFT.  Using  the  three-factor  expression  (3.4.4)  as 
a model  the  radix-r  transform  (where  N = r ) may  be  written  concisely  as 


G = { ir  (I  x M x I ) (D 


p=k  r^ 


x I _) }. 


N/rp  r1*  , r N/r5 


Vr  9 


(3.4.6) 


Here, 


(i)  The  permutation  matrix  r is  given  by 

A = (P,  x I . )(P  x I . o) (P„ / o x I )P„. 

N,r  l,r  N/r  r,r  N/rz  N/r2,r  r N/r,r 

(3.4.7) 

This  permutation  reorders  the  input  elements  into  the  sequence  defined  by 
reversing  the  digits  of  their  original  sequence  position  expressed  as  a 
number  to  base  r. 


1 
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(ii)  The  second  factor  in  the  continued  product  of  eqn  (3.4.6)  is 

always  diagonal  and  consists  of  the  r**  x r15  matrix  D •.  repeated 

rF  »r 

N/r^  times  along  the  diagonal.  As  before,  the  (n^  + m^  ^N/r^  ^)th 

element  of  D D_i  is  (W  j”*®110  where  . 

N/rP  x,r  rp  • 

n = 0,1,  •••  H/z9'1-! 
o 

♦ 

m =0,1,  •••  r-1 
o 

The  only  operations  involved  in  each  pass  or  stage  of  a radix-r 
FFT  are  the  weighting  of  the  input  vector  to  each  pass  by  a diagonal 
matrix  and  the  calculation  of  N/r  r-point  DFT's.  The  only  differences 
from  pass  to  pass  are  the  values  of  the  weights  and  the  selection  of  the 
r elements  upon  which  to  perform  the  r-point  DFT's. 

The  formal  number  of  operations  required  for  a radix-r  FFT  (N^r  ) is 

Nrk  = Nr  log  N 
r 

The  improvement  factor  over  direct  evaluation  of  the  DFT  is  % 

N2/(Nr  logrN)  = (N/log  N) (log  r/r) 

li 

and  is  thus  proportional  to  log  r/r.  The  values  of  this  function  are  shown 
in  the  following  table  : 


r 

log  r/r 

2 

0.347 

3 

0.366 

4 

0.347 

5 

0.322 

6 

0.299 

7 

0.278 

8 

0.260 

9 

0.244 

10 

0.230 

The  function  has  a maximum  for  r=3  with  the  values  for  r=2  and  r=4 
(which  are  ec;ual)  being  slightly  smaller.  For  increasing  values  of  r 
the  function  decreases  slowly.  This  suggests  that  radix-3  is  formally  ^ 

the  most  efficient.  However,  the  radices  2 and  4 provide  significant 

. 
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simplification  to  the  addressing  when  the  FFT  is  performed  on  a binary 
computer. 


More  importantly,  however,  both  radix-2  and  radix-4  algorithms  may 
be  considerably  simplified  by  taking  advantage  of  the  2 and  4 fold 
symmetries  of  the  sine  and  cosine  functions. 


As  stated  previously  a radix-r  FFT  involves  many  r-point  DFT's  at 
each  pass.  The  advantage  of  radix-2  and  radix-4  algorithms  lies  in  the 
fact  that  both  a 2-point  and  a 4-point  DFT  may  be  performed  using  only 
additions  and  subtractions  - no  multiplications  are  required.  This  can 
be  seen  as  follows  :- 


M2  " 

V 

W2 

W°’ 

W2 

= 

1 

1 

r 

S 

NJ  0 

w1 

2_ 

1 

-1 

since 

w2  = 

exp 

(-2iri/2)  = 

-1 

Furthermore , 


(3.4.8) 


fw° 

W4 

w° 

4 

w° 

W4 

W°1 

4 

= 

1 

1 

1 

1 

W° 

W4 

w1 

"4 

w 2 
W4 

w3 

W4 

1 

-i 

-1 

i 

w° 

W4 

w2 

W4 

w4 

W4 

w6 

W4 

1 

-1 

1 

-1 

l 

S 

o 

w3 

W4 

w6 

W4 

w9 

4. 

1 

i 

-1 

-i 

since  = exp  (-2ni/4)  =-i 


As  no  multiplications  are  involved  in  the  r-point  DFT's,  the  only 
multiplications  required  are  the  initial  weighting  of  the  input  vector  to 
each  pass  by  the  diagonal  matrix.  These  weightings  or  phase  changes  to  the 
input  elements  to  each  pass  have  been  called  "twiddling"  the  data.  Consequent- 
ly the  radix-2  and  radix-4  algorithms  are  known  as  the  twiddle  factor 
algorithms. 


F/G  12/1 
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As  an  example,  an  8-point  radix-2  FFT  is 


G - Ma  g 


<I4  X M2)  (°4,2  X V* 


Pass  3 


(I2  x M2  x I2)(D2,2  X I2)*  PaSS  2 


(3.4.10) 


(M2  x I4) (D1>2  x I4) . 


*8,2  9 


Pass  1 


Bit  Reverse 


The  matrices  involved  in  this  expression  sure  shown  in  full  in  Fig. 3.6. 

There  is  a 1:1  correspondence  between  the  matrix  representation 
and  the  commonly  used  signal  flow  graph  of  the  same  transform.  The  signal 
flow  graph  for  the  transform  (3.4.10)  is  shown  in  Fig. 3. 7.  It  should  be 
noted  that  in  all  matrix  expressions  the  "signal  flow"  is  from  right  to 
left.  Each  successive  matrix  multiplication  of  Fig. 3.6  starting  from  the 
right  yields  a vector  of  intermediate  results.  The  elements  of  this 
intermediate  vector  correspond  to  each  successive  vertical  line  of  nodes 
in  the  signal  flow  graph  starting  from  the  left.  If  the  original  data 
are  supplied  in  normal  order  the  operation  of  rearranging  the  elements  into 
bit-reversed  order  is  given  from  (3.4.7)  as  the  product  of  permutation 
matrices  : 


Ag  j = (P-l  2 x I4)  (P2,2  X V (P4f2  X Il) 


(3.4.11) 


It  might  appear  from  the  preceding  discussion  that  it  is  irrelevant 
whether  the  radix-2  or  radix-4  algorithm  is  used.  However,  the  twiddle 
factor  algorithms  involve  N multiplications  in  each  of  the  k passes  and, 
since  there  are  half  as  many  passes  in  radix-4  as  in  radix-2,  there  is  an 
obvious  advantage  in  using  the  radix  4 algorithm.  When  N is  an  odd  power 
of  2 it  is  most  efficient  to  perform  as  many  radix-4  passes  as  possible 
and  one  radix-2  pass.  These  algorithms  are  known  as  "radix-4  + 2". 
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INPUT  ELEMENTS 


a)  Vector  Grouping 


Matrix  Grouping 


go  g3  g6  g9 


gl  g4  97  g10 


g2  95  98  911 


r2  rows 


index  m 


r^  columns 


index  m. 


OUTPUT  ELEMENTS 


go 

4 groups  of  3 

r 

G 

o 

gl 

group  index  m^ 

G1 

g2 

index  within  group  m^ 

G2 

g3 

m = m + 3m, 
o 1 

G3 

g4 

m =0,1,2 
o 

°4 

g5 

= 0 f 1 1 2 , 3 

G5 

g6 

G6 

g7 

G7 

g8 

G8 

g9 

G9 

gio 

G10 

gll 

G11 

_ 1 

- - 

3 groups  of  4 
group  index  n 


index  within  group  n 


n = n + 4n, 
o 1 


n = 0,1, 2, 3 
o 


n1  = 0,1,2 


G G,  G„ 

o 4 8 


G1  G5  G9 


G2  G6  G10 


G3  G7  G11 


r^  rows 


index  n 


r2  columns 


index  n. 


FIG.  3.1  DATA  INDEXING  FOR  FFT  s N=12=4X3 


FIG.  3.7  SIGNAL  FLOW  GRAPH  OF  AN  8-POINT  RADIX  2 COOLEY-TUKEY  FFT 
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CHAPTER  4 : TIME  SERIES  ANALYSIS 

4.1  Introduction 

A "time  series"  in  its  broadest  sense  is  simply  any  function  of 
time,  but,  in  practi-  i,  the  term  is  usually  applied  to  a function  whose 
value  at  any  time  is,  at  least  in  part,  random.  In  part  the  theory  has 
developed  out  of  attempts  to  extend  procedures  valid  for  deterministic 
functions  to  random  functions.  This  generally  means  that  the  methods 
developed  can  legitimately  be  applied,  as  a special  case,  to  purely 
deterministic  processes.  A common  case  of  interest  is  when  the  function 
consists  of  a deterministic  part  combined  with  a random  part  ("noise"). 

Generally  the  aim  of  time  series  analysis  is  to  extract  the 
maximum  information  about  the  deterministic  part  and/or  to  find  the 
statistical  characteristics  of  the  random  part.  The  tools  of  time  series 
analysis  involve  a combination  of  Fourier  methods  and  the  methods  of 
probability  theory.  In  the  following,  several  results  from  probability 
theory  will  be  used  with  little  conment;  these  can  all  be  found  in 
standard  texts  such  as  those  listed  at  the  end  of  this  chapter,  refs.  (1) , 
(2)  and  (3).  Also  included  in  the  references  are  several  works  on  Time 
Series. 

Consider  an  experiment  in  which  some  quantity  g(t)  is  measured. 

The  particular  time  history  observed,  g^t)  say,  is  called  a realisation 
of  the  process.  If  the  experiment  is  repeated  under  identical  conditions 
then  a second  realisation  g2(t)  will  be  observed.  Further  experiments 
give  rise  to  further  realisations  and  the  whole  set  of  realisations  is 
termed  an  ensemble.  In  a non-rigorous  sense  we  may  say  that  if  all  the 
realisations  are  identical,  then  g(t)  is  a deterministic  process,  whilst 
if,  for  each  t,  g(t)  is  a random  variable,  then  g(t)  is  a random  or 
stochastic  process.  (A  rigorous  definition  must  consider  the  cases  when 
(a)  the  g. (t)  are  random  variables  for  some  values  of  t and  identical 
for  others,  and  (b)  the  realisations  are  all  identical  except  for  a few 
exceptional  realisations  that  may  occur  with  zero  probability. 

It  is  possible  to  treat  all  processes  as  stochastic  processes,  the 
deterministic  ones  simply  being  those  in  which,  for  each  time  t,  there  is 
zero  scatter  in  the  values  of  the  time  series.  This  would,  however,  be  a 
rather  wide  definition  of  a deterministic  process.)  The  expected  value 


I 


of  some  function  of  g(t)  is  simply  the  average  over  all  realisations 
(the  ensemble  average)  of  that  function.  Since  for  each  time,  t,  g(t) 
is  a random  variable,  the  process  defines  a set  of  finite  dimensional 
probability  distributions 


1-D 

Ftx^-t^ 

2-D 

F(xirx2;t 

3-D 

F(x  ,x2,x. 

etc 

3-D  F(x1,x2,x3,*t1,t2,t3)  = p[g(t1)<  xx,  g(t2)<  x2,g(t3)<  x^ 


These  finite  dimensional  distributions  are  also  sufficient  to  define  the 
random  process,  according  to  a theorem  of  Kolmogorov. 

A stationary  stochastic  process  is  one  whose  probability  structure 
does  not  change  with  time.  Thus  a strictly  stationary  process  will  be  one 
for  which 

F(x1fx2,  •••;  t^,t2,»*»)  = Ftx^x^,  •••;  t^+T , t2+T  ) 

is  true  for  all  values  of  t.  (In  other  words  the  distributions  depend 
only  one  the  time  differences  t2"ti'  t3_ti  etc-)  However  there  are  some 
parts  of  the  theory  which  can  be  developed  from  a less  restrictive 
assumption.  A weakly  stationary  process  of  order  l is  one  for  which  the 
multivariate  moments 

l i.  I 

e (gft^  1 g(t2)  2-**  g(tR)  n) 

uo  to  order  £ = £.  + £„+•••+£  depend  only  on  the  time  differences 

1 2 n 

t - t^  etc.  A weakly  stationary  process  of  order  2 is  the  most  commonly 
considered  case.  If  the  probability  distribution  associated  with  any  set 
of  times  is  a multivariate  normal  distribution,  the  process  is  called  a 
Gaussian  process.  Since  a normal  distribution  is  fully  characterised  by 
moments  of  first  and  second  order,  in  this  case  second  order  stationarity 
is  enough  to  produce  strict  stationarity. 


V.W-H  «***■ 
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An  ergodic  process  is  one  in  which  all  the  different  realisations 
of  a process  are  adequately  represented  by  different  parts  of  any  one 
realisation.  In  other  words  an  ensemble  average  of  any  function  of  g(t) 
is  equal  to  the  time  average  of  that  function. 


T-*» 


T 

E{<t>Fg(t)J  } = Lim  (1/2T)  / <fr[g(t+t1)]dt1 
-T 


An  ergodic  process  is  necessarily  stationary  and  values  of  the  process  which 
cure  sufficiently  separated  in  time  must  be  uncorrelated. 


4.2  Results  from  Fourier  Transform  Theory 


For  convenience  the  main  results  from  Fourier  transform  theory 
required  in  the  following  have  been  summarised  below: 


a)  Fourier  In  a l_Tr an s fo rm 


Ti 


Time  series 

9 

g(t) 

' C,\<JM\2  dt  < * 

(4.2.1) 

% 

Transform 

9 

G(f) 

= J m g (t)  exp  (-27rif  t)dt 

(4.2.2) 

Inverse 

9 

. g(t) 

= /ooG(f)erp(2irift)df 

(4.2.3) 

Orthogonality 

9 

/ exp(2irift)dt  = 6(f) 

—00 

(4.2.4) 

Convolution 

; 

9x(t) 

* g2  (t)  “/"g^s)  g2(t-s)ds 

(4.2.5) 

b)  Discrete^ 

Fourier  Transform  (N  points) 

Time  series 

l 

g(t) 

: g(t+N)  = g(t) 

(4.2.6) 

Transform 

9 

G(n) 

Hr1 

= (1/N)  I g(t)exp(-2irint/N) 

(4.2.7) 

| 

t=o 

N-l 

Inverse 

1 

g(t) 

= I G(n)exp  (2-rrint/N) 
n=o 

(4.2.8) 

* 

N-l  1 if 

r=n  modulo  N 

Orthogonality 

} 

(1/N) 

E [exp  2nit(r-n)/ig  = 
t=o  0 otherwise 

♦ 

r 


(4.2.9) 


I 
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N-l 


Convolution  ; g^(t)*g2(t)  = (1/N)  E g^(s)  g 2 (t-s)  (4.2.10) 


s=o 


(A  somewhat  abbreviated  but  obvious  notation  has  been  used  here.) 


All  the  above  results  have  been  derived  previously.  One  further 
result  which  can  be  called  the  "multiplication  theorem"  is  needed.  This 
states  that 


(a)  for  the  Fourier  Integral  case 


/“  gL  (t)  g2  (t)  dt  = fjzl  (-f)  G,  (f)  df 


(4.2.11) 


and 


(b)  for  the  DFT  case 


(1/N)  E g1  (t)  g2  (t)  = E G^-n)  G2  (n) 


(4.2.12) 


The  proof  for  the  DFT  case  follows  (that  for  the  integral  being  parallel) 


(1/N)  E g1  (t)  g2  (t) 


(1/N)  E g (t)  {EG  (n)  exp  (2irint/N) } 
t 1 n 


(1/N)  E { E g (t)  exp  (2irint/N) } G2  (n) 
n t 1 


= EG.  (-n)  G_  (n) 
n 1 2 


(4.2.13) 


4.3  The  Power  Spectrum 


4.3.1  An  Invariance  Property  of  Fourier  Transforms 


Consider  the  effect  of  a time  shift  on  the  Fourier 
transform  of  g(t). 


I 


The  Fourier  integral  transform  of  g(t-a)  is 


G"  (f)  - fCn g(t-a)  exp£-2irif  (t-a)J  exp(-2xifa)  d(t-a) 


= / g(t)  exp  (-2xift)  exptfiirifa)  dt 


G (f ) expf€irifa) 


(4.3.1) 


The  two  quantities  G(f),  G' (f ) have  the  same  modulus  but  different  phases. 
Alternatively  it  may  be  said  that  the  quantity 


G*  (f ) G (f ) = | G (f ) |” 


(4.3.2) 


is  invariant  under  shifts  in  the  time  scale.  (The  asterisk  denotes  the 
complex  conjugate.)  The  quantity  P(f)  is  called  the  power  spectral 
density  and,  because  of  the  above  invariance  property,  is  a key  tool  in 
time  series  analysis. 


It  might  be  noted  that  when  g is  real,  P is  an  even  function 
of  f.  In  the  past  it  was  generally  analysed  using  a DFCT. 


4.3.2  Interpretation  of  Power  Spectral  Density 


Suppose  in  the  multiplication  theorem  for  the  DFT,  as  given 
by  (4.2.12),  one  sets 


g2(t)  = g(t) 


g (t)  = g* (t)  so  Gj^-n)  = G*(n) 


then 

(1/N)  IgMt)  g(t)  = (1/N)  l|g(t)|2  = EG*  (n)  G(n)  = EP(n)  (4.3.3) 


This,  of  course,  is  just  Parseval's  theorem.  However  if  g(t)  is 
interpreted  as  the  current  through  a unit  resistance,  the  mean  square  of 
g(t)  is  the  power  dissipation.  This  is  the  reason  for  the  term  power 
spectral  density,  for  P (n)  represents  that  portion  of  the  power  which  is 
contributed  by  the  n'th  frequency  band.  Thus  a signal  having  a periodic 
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component  of  frequency  n/(NAt),  where  At  is  the  sampling  interval, 
will  have  a spike  in  the  spectrum  at  n. 

4.4  Auto  Covariance  and  Auto  Correlation 

4.4.1  Definition 

One  way  of  looking  for  periodic  components  in  a signal  is 
via  the  "auto-covariance".  If  a component  of  period  k At  is  present  in 
the  signal  then  there  will  be  a degree  of  correlation  between  g(t)  and 
g ( t+k) . This  idea  leads  to  the  definition  of  the  auto-covariance  by 

a)  Fourier  integral  case 

C (t)  = /”  g*(s)  g(s+t)  ds  (4.4.1) 

b)  DFT  Case 

N-l 

C (t)  = (1/N)  E g* (s)  g(s+t)  (4.4.2) 

s=o 

For  many  purposes  it  is  useful  to  define  an  auto-correlation  function 
R (t)  by 

R(t)  = C(t)  / C(o)  (4.4.3) 

Note  that 

(i)  C(o)  is  just  the  mean  square  value  of  g(t)  and  so 

C (o)  =£p(n)  (4.4.4) 

(ii)  C(-t)  = C*  (t)  (4.4.5) 

so  when  g and  hence  C is  real,  C is  an  even  function  of  t. 

For  some  purposes  it  is  also  convenient  to  define  a structure 
function  D(t)  given  by 


a) 


Fourier  Integral  Case 


D(t)  = /"|g(s+t)  - g <s) | 2 ds  (4.4.6) 

b)  DFT  Case 

N-l  _ 

D(t)  = (1/N)  £ |g(s+t)  - g(s) | (4.4.7) 

3=0 

It  can  be  shown  that 

D(t)  = 2 C(o)  - [C(t)  + C*(t)j  (4.4.8) 

4.4.2  Wiener  Khintchine  Theorem 

If  the  power  spectral  density  of  a stationary  random 
process  exists,  then  it  is  the  Fourier  transform  of  the  auto-covariance. 

We  will  not  prove  this  theorem  for  the  general  case  of  a stationary  random 
process,  but  only  for  the  particular  DFT  case  considered  previously. 

By  definition  the  transform  of  the  auto  covariance  is 

(1/N)  £ C(t)  exp(-2irint/N) 
t 

- (1/N)  £ {(1/N)  £ g*(s)  g(s+t)}  exp (-2irint/N) 
t s 

I g*(s)  exp (2irisn/N) } {(1/N)  £ g(s+t) 

s t 

exp  [-2iri  (s+t)  n/Nj  } 

= G*(n)  G(n) 

= P(n)  (4.4.9) 

This  result  can  also  be  established  using  the  multiplication  theorem  and 
taking 

gx(t)  - g(t) 

g2(t)  = g* (t-k)  so  G2  (-n)=GA(n)  exp(2rikn/N) 
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Dies  of  Time  Series 


4.5.1  Deterministic  Cases 


Sine  Have,  frequency  f. 

The  spectrum  consists  of  two  lines  at  frequencies  ±f. 


Two  or  more  sine  waves,  different  frequencies.  The 


spectrum  consists  of  discrete  lines  at  ±f^,  ±f2  etc. 


Note;  If  two  frequencies  are  incommensurable 
magnitudes,  (e.g.  f,  /2f)  then  the  time  series  is  not 
periodic.  Nevertheless  the  spectrum  is  still  a line 
spectrum  and  the  time  series  is  called  "quasi-periodic" 


Damped  sine  wave. 

Greater  damping  gives  broader  spectral  peak. 


An  impulse. 

The  power  spectrum  of  a Dirac  Delta  function  is  a 
constant  over  all  frequencies. 


4.5.2  Stochastic  Case  1 - Uncorrelated  Random  Variables 


Consider  a (discrete)  time  series 


{g(t) } - {e (o) , e (1) , e(2),  e(N-l)} 


(4.5.1) 


where  the  c(t)  are  independent  identically  distributed  (i  i d)  random 

. . _ ...  j _ 2 • , j i a m a.  i. . ■ . 


variables,  with  zero  mean,  variance  a and  finite  fourth  moment 
i.e.  for  each  t and  s and  using  E to  denote  expectation 


E(e(t)}  = 0 


(4.5.2) 


E{e(t)e(s)}  * o 6(t-s) 


(4.5.3) 


E{e(t)  } - U. 


(4.5.4) 


Taking  the  DFT  of  g(t)  gives 


G(n)  - (J./N)  £ e (t)  exp(-2wint/N) 


(4.5.5) 
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i.  v.  . ..  '.'j'wjjwv  ...w^i<i..iii.i  . i .■  < .? 


so  that 


P(n)  = G* (n)  G(n) 

= (1/N2)  Z Z e*(s)e{t)exp£-2irin(t-s)/N]  (4.5.6) 

t s 

and 

E{P (n)  } = (1/N2)  ZEE  (e*(s)  e (t)Jexp[-2Trin(t-s)/N] 

t s 

= o2/N  (4.5.7) 


This  is  a constant,  independent  of  n.  Such  a stochastic  time  series, 
whose  power  spectrum  is  constant  for  all  frequencies  is  termed  white 
noise.  (Note  that  the  power  spectrum  of  a delta  function  is  also  a 
constant.)  It  can  also  be  shown  that  if  m ^ n : 

e{  [p  (m)  - e(p  (m)  }|  £p  (n) -E{P  (n)  }J  } “0  (4.5.8) 


and 


E{[.P(n)  -E{P(n)}J2}  s (cV)  [l+6  (n)] 

= E{P(n)2}  jl+6(n)]  (4.5.9) 

This  means  that  for  any  two  given  distinct  frequencies  the  estimates  of 
the  power  spectral  density  at  those  two  frequencies  are  a pair  of 
uncorrelated  random  variables,  and  at  any  single  frequency  the  non- 
negative estimates  of  the  power  spectral  density  will  be  distributed  with 
a variance  equal  to  the  square  of  the  mean  (except  at  zero  frequency  where 
the  variance  is  twice  this  value) . 

The  reliability  of  an  estimate  is  usually  indicated  by  the 
"coefficient  of  variation",  C.V.,  where 

C.V.  = standard  deviation  / mean  (4.5.10) 

In  this  case  the  C.V.  is  unity. 


* 


» 
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Moving  Average  (m.a.) 


where  the  e (t)  are  defined  as  in  the  previous  example,  and  the  a(t)  are 
constants.  If  we  suppose  m < N-l  and  that  a(j)  = o for  j > m we  may 
write  this  as  a convolution 


of  the  two  sequences  N a (t) , e (t) , so  that  the  Fourier  transform  is 
given  by 

r "-1  1 

G(n)  = j(l/N)  E N o (t)  exp  (-2irint/N)J  • 


N-l 

1/N)  E e (t)exp  (-2irint/N) 
t=o 


Writing  e = exp(-2flin/N)  it  can  be  seen  that  the  first  bracket  is 
simply  a polynomial 


exp  -2nin(t-s)/N 


and  the  mean  value  of  P is 


In  this  case  the  power  spectrum  is  no  longer  a constant  independent  of 
frequency.  It  has  a shape  which  depends  on  the  polynomial  A.  It  may 
be  shown  that  any  well  behaved  power  spectrum  can  be  approximated  to 
anv  prescribed  degree  of  accuracy  by  a function  of  the  form  A,  and  the 


resultant  values  of  a(t)  give  a moving  average  which  can  be  used  to 
simulate  a time  series.  If  the  time  series  is  known  up  to  time  t , 
then  it  is  also  possible  to  estimate  the  tt t)  for  t < tx  and  use 
these  as  starting  values  to  predict  the  future  time  history.  The 
disadvantage  of  a moving  average  representation  is  that  it  may  require 
a considerable  number  of  terms  (m  large)  to  achieve  a good  fit  to  an 
arbitrary  power  spectrum. 

It  may  be  shown  that  E (P(m)  P (n) } has  the  same  form 
as  in  the  previous  example,  apart  from  the  multiplying  polynomials,  so 
that  once  again  spectral  estimates  at  different  frequencies  are  nearly 
independent,  and  the  coefficient  of  variation  is  again  unity.  Mote  that 
a necessary  condition  for  the  m.a.  process  to  be  invertible  (i.e.  useful 
in  forecasting)  is  that  the  polynomial  A(e)  have  all  its  roots  outside 
the  unit  circle. 

4.5.4  Stochastic  Case  3 - Autoregression  (a.r.) 

An  autoregression  is  a process  where  the  value  at  one 
instant  depends  on  the  values  at  previous  instants  e.g. 

g(t)  =0.5  g(t-l)  +0.1  g(t-2)  + e(t) 

In  general  we  write 


m 

E 6(j)  g (t-j)  = e(t)  (4.5.16) 

j=o 


(Note:  A Markov  process  is  a special  case  of  an  a.r.  for  which  m=l) 

By  defining  8(j)  = 0 for  j > m we  may  again  write  this  as  a convolution, 
and  upon  Fourier  transformation  we  obtain 


jjl/N) 


G(n)  = 


N-l 

E NB(t)exp(- 
t=o 

N-l 

= (1/N)  I e (t)exp(-2irint/N) 
t=o 

[(1/N)  E e(t)exp(-2uint/N)]  /B{«) 


(4.5.17) 


— 
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4.6  Reliability  of  Spectral  Density  Estimates 

If  a signal  is  deterministic  then  one  example  (realisation) of  the 
signal  is  enough  to  define  the  power  spectrum.  * 

On  the  other  hand  the  preceding  three  (very  general)  examples  of 
stochastic  processes  all  have  a coefficient  of  variation  (of  the  spectral  • 

density  estimate  at  any  given  frequency)  of  unity,  which  implies  a very 
high  scatter. 

There  are  basically  two  different  methods  to  improve  the  reali- 
ability  of  spectral  density  estimates. 

a)  Averaging  over_realAs^at.ions 

This  is  the  most  computational ly  economical  method. 

Determine  the  power  spectra  of  several  different  realisations, 
gx(t),  92(t),  etc.,  of  the  process,  and  average  the  resultant 
spectral  density  estimates,  P^n) , P2(n)  etc.  at  each  frequency. 

b)  Averaging  over_freguency 

This  is  a means  of  trading  off  frequency  resolution  for 
improved  realiability  of  estimates.  Since  P(m)  and  P(n)  are 
nearly  independent  variables  if  m / n,  we  may  average  adjacent 
estimates.  (This  of  course  will  only  work  if  E {P (n)}  is 
reasonably  constant  over  frequency  intervals  of  size  m - n.) 

One  procedure  which  comes  into  this  category  is  termed  Hanning. 

Using  a prime  to  denote  the  more  reliable  estimate 

] 

P1  (n)  * h P(n-l)  + P(nj  + h p(n+l) 

This  procedure  frequently  introduces  a much  greater  reliability 

improvement  than  would  be  expected,  for  quite  a small  trade-off 

in  resolution.  This  is  because  a time  series  selected  at  random 

is  likely  to  end  up  at  a value  significantly  different  from  its 

starting  value.  The  assumption  of  periodicity  in  the  data  means 

that  this  abrupt  signal  jump  gives  rise  to  high  frequency 

fluctuations  in  the  spectral  density  so  that  P(n)  and  P(n-l) 

are  negatively  correlated  for  a given  realisation.  Such  negative 

correlations  did  not  show  up  in  the  preceding  exanples  because  the  • 

expected  values  were  averages  over  all  realisations. 
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With  either  method  of  improving  reliability,  the  averaging  of 
"k"  independent  variables  reduces  the  scatter  as  l//k.  To  obtain  95% 
confidence  bands  (approx.  2 standard  deviations)  within  10%  of  the  mean 
(the  standard  figure  for  unthinking  engineers  to  pull  out  of  the  hat) 
will  require  the  averaging  of  400  independent  estimates! 

4.7  Linear  Systems 

The  introductory  part  of  this  section  is  taken  in  abbreviated 
form  from  Solodovnikov' s book  (see  ref. 10)  where  further  details  may  be 
obtained . 

Consider  a linear  system  with  k degrees  of  freedom.  Here  "k" 
independent  co-ordinates  cure  needed  to  specify  the  state  of  the  system 
at  any  instant,  and  "k"  independent  excitations  may  be  applied  to  the 
system.  However  because  of  the  principle  of  superposition  we  may  consider 
only  one  of  these  excitations  at  a time  and  set  all  other  excitations  to 
zero.  Say  the  one  non-zero  excitation  is  g^t),  and  one  of  the  co-ordinates 
giving  the  resultant  state  of  the  system  is  g 2 (t) . The  two  functions  g^  and 
g2  are  related  by  a linear  differential  equation 

(bqDq  + b^_LDq-1  + ...  + bQ)  g2  = (arDr  + a^  Dr_1  + . . . + aQ)  gx 

(4.7.1) 

where  D = d/dt  and  (for  a second  order  system)  q £ 2k,  r ^ 2 (k-1) . 

Bearing  in  mind  that  the  Fourier  transform  of  dg/dt  is  2irifG(f),  it 
follows  that  on  taking  the  Fourier  transform  of  both  sides  of  equation 
(4.7.1),  and  ignoring  the  initial  conditions  of  the  system  (which  is 
equivalent  to  ignoring  the  transient  oscillations)  the  result  is 

(b  eq  + ...  + b ) G_  (f)  — (a  * + • • • + a ) G.  (f)  (4.7.2) 

q o 2 r o i 

where 

e = 2irif 

The  transfer  function  of  the  system  is  given  by 

H (f ) = G2  (f)  / Gx  (f)  (4.7.3) 


A (■)  / B(e) 


(4.7.4) 


where  the  polynomials  A and  B sure  given  by 


*■  t 

A(a)  = £ a s 

t=o 

« t 

B(a)  = £ bt  ® 
t=o 


The  transfer  function,  H(f),  may  be  calculated  if  the  parameters  of  the 
system  are  known,  or  experimentally  observed  by  measuring  the  output,  g^, 
caused  by  a known  input  g^  The  transfer  function  is  usually  presented 
as  two  graphs  showing  its  modulus  and  phase. 

Examples 


(a)  If  the  input  g^t),  is  white  noise,  then  G^f)  is  a constant, 
and  the  Fourier  transform  of  the  output  is  (apart  from  a constant 
multiplier)  equal  to  the  transfer  function. 

(b)  If  the  input,  g^t),  is  an  impulse  (Dirac  Delta),  then  the  output 
g^ (t) , is  termed  the  impulse  response . The  Fourier  transf orm  of  an 
impulse  is  also  a constant,  so  again  the  Fourier  transform  of  the 
output  is  equal  to  the  transfer  function. 

(c)  If  the  input,  g^t),  is  a unit  sine  wave  of  frequency  f,  the 
output  is  also  a sine  wave  of  frequency  f,  with  magnitude  equal 

to  the  amplitude  of  the  transfer  function  at  frequency  f,  and  phase 
(with  respect  to  the  input  wave)  equal  to  the  phase  of  H(f). 


(d)  If  the  input,  g^t),  is  a step  function 


gx(t) 


0 

1 


t < 0 
t > 0 


then  the  resultant  output  is  termed  the  step  response.  The  Dirac 
Delta  is  the  differential  of  the  step  function,  so  the  Fourier 
transform  of  the  Dirac  Delta  (known  to  be  a constant)  is  2irif  times 
the  Fourier  transform  of  a step  function.  Therefore  the  impulse 
response  is  2irif  times  the  step  response. 


As  previously  we  may  show  that 


THE  CROSS  POWER  SPECTRUM  IS  THE  FOURIER  TRANSFORM  OF  THE 


Note 


In  the  case  when  (t)  is  the  input  to  a linear  system,  and 
g (t)  is  the  output  we  have 


Transfer  Function 


g (t)  = exp(2irif  t + i?  ) 


g (t)  = exp(2IIif1t  + iC  ) 


Then 


4.8.2  Coherence  Function 


Consider  two  time  series  (t) , g 2 (t) , for  which  we 
postulate  g2(t)  is  made  up  of  a linear  function  of  g^ (t)  and  some  other 
independent  signal 

g2(t)  = + g2e(t) 

Thus 

G2(f)  = G2L(f)  + G2e(f)  = H(f)  Gl(f)  + G2e(f) 

This  is  now  a simple  linear  relationship,  so  to  determine  the  relative 
importance  of  the  two  (randomly  varying)  components  we  measure  the 
correlation  between  G^f)  and  G2(f).  This  is  called  the  coherence 
function. 

y (f ) = e{|g*  (f)  g2 (f ) | > //[e  {|g1|2}  -e{ |g2 | 2>J 

= E |P12(f)  |//[E{P11(f)}  • E{P22(f)}] 

At  frequencies  where  G^  is  much  greater  than  G-,e, 
converse  is  true  y 0 


y -*■  1;  when  the 
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The  purpose  of  this  chapter  is  to  consider  some  of  the  more 
practical  aspects  of  the  operation  and  use  of  Fourier  transforms. 

The  discussion  will  be  limited  to  the  case  of  equi-spaced  data 
analysed  by  the  Fast  Fourier  transform  algorithm  developed  by  Cooley 
and  Tukey.  When  the  term  "Fourier  transform"  is  used  the  discrete 
Fourier  transform  is  implied.  Before  commencing  the  discussion  on 
actual  computer  techniques  some  of  the  practical  consequences  of  the 
concepts  introduced  in  previous  chapters  will  be  discussed. 

5. 2 Frequency  Resolution  and  Range 

Suppose  that  an  experiment  is  about  to  be  conducted  where  data 
from  one  or  more  parameters  are  to  be  collected  over  a period  of  time 
(or,  possibly,  space)  and  that  the  data  are  to  be  analysed  using  Fourier 
techniques.  After  ensuring  that  a transducer  is  available  to  record  the 
parameters  to  the  desired  accuracy,  the  next  consideration  is  that  of 
the  recording  system.  The  cost  and  complexity  of  the  recording  system 
depend  on  several  factors  including  the  sampling  interval  between  each 
data  value  and  the  length  of  each  data  sampling  run.  These  two  factors, 
in  turn,  depend  on  the  frequency  range  and  resolution  that  is  desired. 

If  N data  values 

{ V V g2  gN-i  } 

are  recorded  at  times 

{ 0,  At,  2At  ...  (N-l)At  } 

i.e.,  at  a sample  interval  of  At,  then  the  N point  Fourier  transform 
of  the  data  may  be  written  as 


{ gq,g1,g2  ...  gn_1  } 


real  part  of  G..  = real  part  of  G^_ ^ 

imaginary  part  of  G.  = - (imaginary  part  of  GNj) 

so  that,  for  all  practical  purposes  the  useful  range  of  the  Fourier 
transform  is  from  Gq  to  GN/2-  Hence,  the  maximum  frequency,  f^. 


able  to  be  determined  is  given  by 


t 
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Applying  this  to  an  example,  assume  that  in  order  to  determine 
a vibration  in  an  aircraft  wing  it  is  necessary  to  obtain  the  Fourier 
transform  of  the  wing  tip  acceleration  for  frequencies  between  0.1  and 
20  Hz  in  0.05  Hz  steps.  Here 


f = 20  Hz  and  Af  = 0.05  Hz 
max 

Hence , 


At  = l/fR  = 1/40  second 

and 

N = f /Af  = 800 

K 


Thus,  in  this  example,  a sample  rate  of  at  least  40  times  per  second 
and  a record  length  of  at  least  20  seconds  is  required.  These  figures 
cure  minimum  requirements  and  may  have  to  be  modified  in  the  light  of 
the  discussion  in  the  next  sections. 


5.3 


The  maximum  frequency  estimate  that  cam  be  computed  using  a DFT 
is  equal  to  half  the  sampling  frequency.  That  frequency  which  is  half 
the  sampling  frequency  is  generally  known  as  the  Nyquist  frequency. 
However  the  data  being  recorded  may  have  oscillations  of  a higher 
frequency  so  it  is  necessaury  to  investigate  the  effects  of  these 
oscillations  on  the  Fourier  transform. 


In  Figs.  5.1(a),  (b)  amd  (c)  are  shown  the  same  sinusoidal 
oscillation  which  has  been  sampled  at  three  different  sample  rates.  In 
the  first  case  there  is  no  ambiguity  in  the  curve  fitted  to  these  data 
using  Fourier  techniques.  However,  in  the  second  amd  third  cases  there 
are  ambiguities  and  the  computed  transform  will  be  that  of  the  oscillation 
which  has  a frequency  between  0 and  the  Nyquist  frequency  (shown  as  a 
dotted  line  on  the  Figure) . 
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In  general  the  DFT  estimate  of  the  true  Fourier  transform  at 
the  frequency  f,  where 

0 < f 4 fR/2 

is  dependent  on  the  value  of  the  true  transform  at  frequencies 

f,  f - f,  f + f,  2fD  - f,  2f  + f,  etc. 

R R R R 

Therefore,  to  get  am  accurate  estimate  of  the  true  tramsform  at  frequency 
f the  amplitude  of  all  oscillations  at  frequencies  greater  than  the 
Nyquist  frequency  must  be  very  much  smaller  than  those  at  frequencies  less 
tham  the  Nyquist  frequency.  Steps  must  be  taken  to  ensure  that  these 
higher  frequencies  are  dealt  with,  either  by  increasing  the  sampling 
frequency  or  by  filtering  the  data  before  sampling  it. 

Sometimes  it  is  possible  to  live  with  aliasing  if  only  one  or  two 
distinct  known  frequencies  appear  above  the  Nyquist  frequency.  However 
these  frequencies  will  be  "folded  back"  over  other  lower  frequencies  and 
may  mask  them  out.  Consider  a 50  Hz  noise  sampled  at  80  Hz.  The  50  Hz 
noise  will  be  folded  back  and  will  appear  at  30  Hz  (its  alias)  and  this 
may  cover  another  as  yet  undetected  oscillation  near  30  Hz. 

The  best  method  of  overcoming  the  problem  of  aliasing  is  to  filter 
the  incoming  data  before  recording  it.  It  is  not  possible  to  design  a 
filter  with  a flat  response  up  to  a cut-off  frequency  and  then  zero  response 
above  it;  the  complexity,  size  and  cost  of  filters  usually  depend  on  the 
sharpness  of  the  cut-off.  In  Fig.  5.2  is  shown  the  response  of  some  of  the 
more  well-known  filters.  Note  that  the  filters  may  change  the  phasing  of 
the  signal  and  that  for  frequencies  just  above  the  cut-off  frequency  the 
energy  may  not  be  reduced  by  a large  enough  amount  to  stop  aliasing.  If 
the  area  of  interest  includes  the  frequency  range  near  the  cut-off 
frequency  it  may  be  necessary  to  correct  the  calculated  Fourier  transform 
using  the  filter's  response. 


A quick  way  to  determine  whether  aliasing  is  present  in  the  data  is 
to  make  a special  record  of  data  at  a higher  sampling  frequency  than  that 
used  initially.  If  the  transform  in  the  region  of  interest  is  modified  as 
compared  with  that  for  the  initial  sampling  then  aliasing  is  present. 


♦ 
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of  Fourier  Estimates 


In  Fig.  5.3(b)  is  shown  the  25*  point  DFT  of  the  data  shown  in 
Fig.  5.3(a).  A glance  at  these  results  shows  that  it  is  rather  difficult 
to  determine  which  of  the  many  peaks  are  meaningful.  It  is  obvious  that 
to  get  any  useful  information  out  of  the  data  some  averaging  or  smoothing 
techniques  will  have  to  be  considered. 

One  method  of  obtaining  more  reliable  estimates  is  to  average  the 
Fourier  estimates  over  a small  frequency  band.  This  will,  of  course, 
reduce  the  frequency  resolution  but  adds  to  the  stability.  Using  a 
weighted  average  such  as  Hanning  i . e . , 


Gj  ’ (-Gj-i * 2 Bj  - Vi> /4 


(5.4.1) 


has  another  beneficial  effect  in  that  the  Hanning  window  reduces  the 
leakage  of  high  energy  caused  by  discontinuities  at  the  end  of  the  data 
records  (See  Section  5.5). 

The  Fourier  transform  is  a complex  function  and,  rather  than 
consider  it  as  a real  and  imaginary  part,  it  is  sometimes  convenient  to 
consider  it  as  a magnitude  and  a phase.  In  many  cases  the  phase  is  not 
required  and  so  the  magnitude  alone  can  be  used.  It  is  usual  to  define 
the  Power  Spectrum 


(5.4.2) 


where  the  asterisk  denotes  a complex  conjugate.  The  power  spectrum  is  a 
real  function  and  allows  the  use  of  other  averaging  techniques. 

If  there  are  more  than  data  values  available  for  analysis  it 

is  possible  to 

(i)  take  blocks  of  Xf>6  points  and  compute  the  DFT  for  each  block, 

(ii)  compute  the  power  spectrum, 

(iii)  average  the  power  spectral  estimate  corresponding  to  each  frequency 


In  Fig.  5.4  is  shown  the  result  when  this  technique  is  applied  to  sake 
20  blocks  of  2£t  points  recorded  during  the  same  data  run  as  the  data 
in  Fig.  5.3(a).  The  difficulties  of  determining  meaningful  information 
from  the  power  spectrum  of  Fig.  5.3(c)  are  absent  in  Fig.  5.4. 

In  general,  if  there  are  N independent  data  values  used  in  the 
calculation  of  m frequency  estimates  of  the  power  spectrum  then  the 
90%  confidence  limits  can  be  calculated  as  (see  ref.  1) 

10  { 2/ (3k  - 1)  + l/v  (k  - 1 )}  db  (5.4.3) 

where  the  number  of  degrees  of  freedom,  k,  is  given  by 

k = 2N/m  (5.4.4) 

Note  that  as  N increases,  k increases  and  the  confidence  bands  reduce 
approximately  as  /n  for  a constant  number  of  frequency  estimates,  so  N 
must  be  very  large  to  get  good  stable  estimates  of  the  power  spectrum. 

In  Figs.  5.3(c)  and  5.4  the  90%  confidence  limits  have  been  shown  by  a 
dotted  line  to  give  some  idea  of  the  effect  of  increasing  the  number  of 
data  values  by  a factor  of  20. 

5.5  Effect  of  Large  Oscillations 

Experimental  data  often  contain  dominant  oscillations  which  may 
swamp  other  smaller  oscillations  of  interest.  In  Section  5.3  the  case 
where  these  oscillations  have  frequencies  larger  than  the  Nyquist 
frequency  was  discussed  and  it  was  concluded  that  these  frequencies  have 
to  be  filtered  out  before  sampling.  However,  the  case  of  large 
oscillations  at  lower  frequencies  can  be  dealt  with  after  the  data  are 
sampled . 

If  there  is  a large  amount  of  energy  in  a small  frequency  band, 
the  Fourier  transform  of  the  data  will  suffer  from  "leakage"  and  a 
diffraction  pattern  effect  may  be  observed . The  side  lobes  of  this 
diffraction  pattern  may  swamp  out  those  details  of  the  transform  which 
may  be  required. 
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It  is  possible  to  overcome  this  problem  by  pre-processing 
the  data  to  reduce  the  contribution  of  the  dominant  frequency  components. 
Then,  after  taking  the  transform,  compensation  for  the  pre-processing 
may  be  made.  For  example,  if  it  is  known  that  there  is  a large  amount 
of  energy  in  the  very  low  frequency  end  of  the  transform,  the  data  could 
be  pre-processed  by  forming  a new  time  series 


{ go'  gr  g2  gN-i 


according  to 


= gj  “ 


(5.5.1) 


It  can  be  shown  that  the  power  spectrum  of  this  new  time  series 


< p;-  p;-  p2  — p»/2  > 


can  be  corrected  using 


P..  = P'  / (1  - cos  (2tt  j/N)  }d 


(5.5.2) 


This  pre-whitening  and  post-darkening  (as  it  is  commonly  known)  will 
reduce  the  effect  of  side  lobes  on  the  high  frequency  side  of  a large 
oscillation  but  will  increase  the  side  lobes  on  the  low  frequency  side. 

In  general,  any  pre-processing  of  the  data  which  produces  an  uncorrected 
power  spectrum  that  is  as  close  as  possible  to  white  noise  (i.e.  a flat 
spectrum)  will  optimally  reduce  the  effect  of  side  lobes  provided,  of 
course,  that  the  function  to  correct  for  the  pre-processing  can  be  found. 

If  it  is  preferred  not  to  use  this  type  of  pre-processing  it  is, 
however,  advisable  at  least  to  subtract  off  the  mean  or  even  the  straight 
line  of  best  fit  ("the  linear  trend")  from  the  data  before  computing  the 
transform.  If  no  post-correction  is  applied,  this  will  have  an  effect  on 
the  first  few  low  frequency  estimates  of  the  transform  but  it  is  usually 
easier  to  live  with  this,  than  to  have  a larger  number  of  estimates 
swamped  by  the  side  lobes  of  very  low  frequency  oscillations.  If  no 
correction  for  low  frequency  oscillations  is  applied  the  transform  can 
easily  be  subject  to  misinterpretation. 
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5. 6 Computational  Requirements 

Using  the  information  given  in  the  previous  sections  it  should  be 
possible  to  come  to  a decision  on  the  sampling  interval,  duration  of 
data  runs,  filtering  and  pre-processing  required.  However,  before  the 
experiment  is  attempted  it  is  useful  to  look  at  the  processing 
requirements  of  the  Fourier  transform  itself. 

For  all  but  the  case  of  very  small  amounts  of  data,  the  Fast  Fourier 
Transform  algorithm  is  much  the  faster  method  of  analysis  so  the  discussion 
will  be  confined  to  it.  The  computer  to  be  used  must  have  the  storage 
requirement  of  the  N data  values  to  be  transformed  plus  extra  storage  for 
the  program  and  data  manipulation.  In  general,  the  more  data  manipulation 
space  available,  the  faster  the  computation  can  be  made. 

A typical  32  point  transform  using  the  Sande  modification  of  the 
Cooley-Tukey  method  is  shown  step  by  step  below.  This  routine  is  written 
to  "inplace"  compute  the  transform,  i.e.  it  puts  each  calculation  back 
into  the  same  data  storage  location  and  then  sorts  them  out  at  the  end  of 
the  calculation.  (As  there  are  many  proven  versions  of  FFT  computer 
programs  already  available,  it  is  generally  advisable  to  use  one  of  the 
existing  programs,  possibly  with  modifications  for  individual  needs, 
rather  than  to  attempt  to  write  one  from  scratch.) 

The  routine  shown  below  starts  with  the  complex  data  to  be  trans- 
formed in  locations  0 to  31.  Each  location  contains  2 numbers  - the  real 
and  imaginary  part  of  the  data.  In  this  case  the  data  to  be  transformed 
is  real. 

The  total  transform  is 
A-l 

G(n  + m c + HBC)  * { Z exp  (2lIial/A) . exp(2Hiam/AB) . 
a=o 

B-l 

{ I exp  (2iribm/B) . exp  (2iri  (a+hA)/ABC) . 
b«o 

C-l 

{ I exp  (2iric/C).  g (a+bA+cAB) }}}  (5.6.1) 

c*o 


where  N = ABC  with  A,  B,  C integers 


♦ 


and  0 $ l < A-l 
0 4 m 4 B— 1 
0 4 n 4 C-l 

Using  a Radix  4 4-2  transform, 

C = 4,  B = 2 and  A = N/8 


and  (5.6.1)  becomes 

N/8-1 

G(n  + 4m  + 81  ) = { I exp(16niaJl/N)  .exp  (8friam/N) . 

a=o 

1 

{ E exp  (2irin  (8a  + bN)/8N). 
b=o 

3 

{ E exp(2IIicn/4) . g(a  + Nb/8  + Nc/4)}}} 

c=5°  (5.6.2) 

Using  the  fact  that 


exp  (2iriX)  = 1 

if 

X 

is 

integral 

- -1 

if 

X 

is 

1/2,  3/2,  5/2  ... 

= +i 

if 

X 

is 

1/4,  5/4,  9/4  ... 

= -i 

if 

X 

is 

3/4,  7/4,  11/4  ... 

this  can  be  further  simplified. 


Defining 

F_j  - cos(2w}/32)  + i sin(2irj/32)  (5.6.3) 

a step-by-step  32  point  transform  is  shown  in  Figs.  5.5  to  5.10.  The  values 
of  F.  used  in  the  calculation  axe  shown  in  the  Table  below. 


The  first  calculation  (Step  1),  which  is  of  radix  4 requires  factors 
derived  from  2n/32;  all  other  steps  only  involve  factors  derived  from 
2ir/16.  Hence,  after  Step  1,  only  the  even  number  factors  in  the  Table 
cure  used. 


j 

cos  (2tt  j/32) 

sin(2irj/32) 

0 

1.0000 

0.0000 

1 

0.9808 

0.1951 

2 

0.9239 

0.3827 

3 

0.8315 

0.5556 

4 

0.7071 

0.7071 

5 

0.5556 

0.8315  T 

6 

0.3827 

0.9239 

7 

0.1951 

0.9808 

8 

0.0000 

1.0000 

10 

-0.3827 

0.9239 

12 

-0.7071 

0.7071 

14 

-0.9239 

0.3827 

Table  : Factors  used  in  32-Point  FFT 


* 


j 

I 

y 

t 


I 


Step  6 (Fig.  5.10)  is  the  most  complicated  and  is  usually  calculated 
with  the  aid  of  a bit  reversal  routine,  sometimes  written  in  machine 
language  to  increase  speed.  The  total  routine  uses  approximately  1.3N 

» 

storage  locations  and  it  uses  the  technique  of  dividing  each  calculation 
by  2 as  the  calculation  proceeds.  Using  this  technique  the  magnitudes 
of  the  numbers  throughout  the  calculation  remain  approximately  the  same,  • 

minimising  round-off  errors  and  allowing  fixed  point  variables  to  be  used 

if  this  is  a constraint  posed  by  the  computer  used  in  the  analysis. 

' 

.1 

1 


( 
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5.7  Calculation  of  Convolutions 


Before  closing,  brief  reference  will  be  made  to  the  calculation 
of  convolutions  using  the  FFT.  Here  it  is  a question  of  evaluating  a 
sum  such  as 

N-l 

z=  l y • h m = 0,1,  ...  N-l  (5.7.1) 

m . 1 m- 1 

j=o 

Reference  to  the  precautions  necessary  to  avoid  wrap-around  error  has 
already  been  made  in  Chapter  2 and  this  will  not  be  repeated  here.  How- 
ever, it  is  worth  remarking  that,  because  of  the  speed  of  the  FFT 
algorithm the  convolution  (5.7.1)  is  usually  calculated  in  an  indirect 
manner  as  follows  : 

(i)  From  the  data  values  y and  h compute  the  transforms 

m m 

Y and  H (using  the  FFT) , 
n n 

(ii)  By  virtue  of  the  convolution  theorem,  the  transform 

Z of  z is  then  obtained  from 
n m 

Z = Y H (5.7.2) 

n n n 

(iii)  Compute  z by  applying  the  inverse  transform  to  eqn 

xn 

(5.7.2). 

This  is  generally  a faster  procedure  than  a direct  evaluation  of 
eqn  (5.7.1). 
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(b)  Fourier  transform  of  data;  real  part  above, 
imaginary  part  below  (function  of  frequency) 


(c)  Power  spectrum  (function  of  frequency) 


FIG.  5.3  FOURIER  ANALYSIS  OF  DATA  WITHOUT  AVERAGING 


» 


FIG.  5.4  FOURIER  ANALYSIS  OF  DATA  WITH  AVERAGING; 
POWER  SPECTRA  OBTAINED 

(a)  From  20  summations 

(b)  From  200  summations 

OF  256  POINT  FOURIER  TRANSFORM 
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* OLD  ARRAY  FACTOR  NEW  ARRAY 


REAL 

I MAG 

REAL 

IMAG 

0 50.00000 

0.00000* 

^50.00000 

0.00000 

1 86.60254 

0.00000  \ 

*43. 30127 

0.00000 

2 100.00000 

0.00000  \ 

/ 

25.00000 

0.00000 

3 86.60254 

0.00000  \ 

/ 

-0.00000 

0.00000 

4 50.00000 

0.00000  \ 

/ 

-25.00000 

0.00000 

5 0.00000 

0.00000  \ 

/ 

-43.30127 

0.00000 

6 -50.00000 

0.00000  \ 

/ 

-50.00000 

0.00000 

7 -86.60254 

0.00000  \ 

/ 

-43.30127 

0.00000 

8-100.00000 

0.00000  { ^ 

>-25. 00001 

0.00000 

9 -86.60254 

0,00000  \ 

V / 

-0.00000 

0.00000 

10  -50.00000 

0.00000  \ J 

A / 

25.00000 

0.00000 

11  0.00000 

0.00000  \ / 

\ / 

43.30127 

0.00000 

12  50.00000 

0.00000  W 

\ / 

50.00000 

0.00000 

13  86.60254 

0.00000  X 

\ / 

43.30127 

0.00000 

14  100.00000 

0.00000  /\ 

Y 

25.00000 

0.00000 

15  86.60254 

0.00000  / \ 

0.00000 

0.00000 

16  50.00001 

0.00000-4 A, 

-/-■ «»F0 

-0.00000 

0.00000 

17  0.00000 

0.00000 

V FI 

42.46924 

8.44765 

18  -5 0.00000 

0.00000  J 

rV  F2 

69,29094 

28.70123 

19  -86.60254 

0.00000  / 

\ F3 

72.00734 

48.11376 

20-1 00.00000 

0.00000  / 

\ F4 

53.03297 

53,03297 

21  -86.60254 

0.00000  / 

\ F5 

24.05688 

36.00367 

22  -50.00000 

0.00000  / 

\ F6 

0.00000 

0.00000 

23  0,00000 

0,00000  / 

fc  F7 

-8.44765 

-42.46924 

24  49.99999 

o.ooooo  L 

♦ FO 

0,00000 

-75.00000 

25  86.60254 

0.00000 

FI 

16.89530 

-84.93847 

26  100.00000 

0.00000 

F2 

28.70123 

-69.29094 

27  86.60255 

0.00000 

F3 

24.05688 

-36.00367 

28  50.00001 

0.00000 

F4 

0.00000 

-0.00000 

29  0.00000 

0.00000 

F5 

-36.00367 

24.05688 

30  -50.00001 

0.00000 

F6 

-69.29094 

28.70124 

31  -86.60254 

0.00000 

F7 

-84.93848 

16.89530 

FOR  J - 0 TO  7 

G(J>=CG(J)+G(J+ 16)3/2 
G ( J+8 ) = CG  ( J+8 ) +G ( J+24 ) 3/2 
G< J+16)=F( J)*CG< J)-G( J+16)3/2 
G < J+24 ) =F ( J ) *Z G ( J+8 ) -G ( J+24 ) 3/2 

STEP  1 OF  A 32  POINT  FOURIER  TRANSFORM. 


FIGURE  5.5 


* OLD  ARRAY  FACTOR  NEW  ARRAY 


REAL 

I MAG 

REAL 

IMAG 

0 

50.00000 

0.00000  ^ — 

3 

r12. 50000 

0.00000 

1 

43.30127 

0.00000  \ 

y 

21.65063 

0.00000 

2 

25.00000 

0.00000  \ 

25.00000 

o.qoooo 

3 

-0.00000 

0.00000 

\ / 

21.65064 

0.00000 

4 

-25.00000 

0.00000 

XX 

12.50000 

0.00000 

5 

-43.30127 

0.00000 

X 

0.00000 

0.00000 

6 

-50.00000 

0.00000 

x\ 

-12.50000 

0.00000 

7 

-43.30127 

0.00000  y 

\ 

-21.65063 

0.00000 

8 

-25.00001 

0.00000 

^-FO 

37.50000 

0.00000 

9 

-0.00000 

0.00000 

F2 

20.00257 

8.28533 

10 

25.00000 

0.00000 

F4 

0.00000 

0.00000 

11 

43.30127 

0.00000 

F6 

-8.28533 

-20.00257 

12 

50.00000 

0.00000 

F8 

0.00000 

-37.50000 

13 

43.30127 

0.00000 

F10 

16.57066 

-40.00514 

14 

25.00000 

0.00000 

F12 

26.51648 

-26.51648 

15 

0.00000 

0.00000 

. F14 

20.00257 

-8.28533 

16 

-0.00000 

0 , 00000  r — 

*7 

-0.00000 

-37.50000 

17 

42.46924 

8.44765  \ 

/ 

29.68227 

-38.24541 

18 

69.29094 

28.70123  \ 

v / 

48.99609 

-20.29485 

19 

72.00734 

48.11376 

\ / 

48,03211 

6.05504 

20 

53.03297 

53.03297 

26.51649 

26.51648 

21 

24.05688 

36.00367 

SX 

-5.97340 

30.03027 

22 

0,00000 

0.00000 

X \ 

-34.64547 

14.35062 

23 

-8.44765 

-42.46924  / 

-46.69306 

-12.78697 

24 

0.00000 

-75 « 00000  — 

— JuFO 

-0.00000 

37.50000 

25 

16.89530 

-84 . 93847 

F2 

-6.05503 

48.03210 

26 

28.70123 

-69.29094 

F4 

-20.29482 

48.99606 

27 

24.05688 

-36.00367 

F6 

-29.68225 

38.24538 

28 

0.00000 

-0.00000 

F8 

-26.51649 

26.51648 

29 

-36.00367 

24.05688 

F10 

-17.01077 

25.45843 

30 

-69.29094 

28.70124 

F12 

-14.35062 

34.64544 

31 

-84.93848 

16.89530 

FI  4 

-23.97523 

42.05870 

FOR  K = 0 , 16 

FOR  J = 0 TO  7 

G< J+K)=rG( J+K)+G< J+K+8) 3/2 
G ( J+K+8 ) =F ( 2* J > *CG  < J+K ) -G ( J+K+8 ) 3/2 

STEP  2 OF  A 32  POINT  FOURIER  TRANSFORM ♦ 


FIGURE  5.6 


* 

OLD  ARRAY 

FACTOR 

NEW 

ARRAY 

REAL 

IMAG 

REAL 

IMAG 

0 

12.50000 

0. 00000  xr 

> 

12.50000 

0.00000 

1 

21.65063 

0.00000  X. 

10.82532 

0.00000 

2 

25.00000 

0.00000  X 

6.25000 

0.00000 

3 

21.65064 

0.00000 

V 

0.00000 

0.00000 

4 

12.50000 

0.00000 

— a^Fo 

-0.00000 

0.00000 

5 

0.00000 

0.00000 

F4 

7.65465 

7.65465 

6 

-12.50000 

0.00000  - 

F8 

0.00000 

18.75000 

7 

-21.65063 

0.00000 

F12- 

15.30930 

15.30930 

8 

37.50000 

0.00000  ^ 

18.75000 

-18.75000 

9 

20.00257 

8.28533  \ 

18.28662 

-15.85990 

10 

0.00000 

0.00000  X 

13.25824 

-13.25824 

11 

-8.28533 

-20.00257 

5.85862 

-14.14395 

12 

0.00000 

-37.50000  

^FO 

18.75000 

18.75000 

13 

16,57066 

-40.00514 

F4  - 

15.85989 

18.28661 

14 

26.51648 

-26.51648 

F8  - 

13.25824 

-13.25824 

15 

20.00257 

-8.28533 

F12 

14.14394 

-5.85861 

16 

-0.00000 

-37. 50000  x- 

-13  ♦ 25824 

-5.49176 

17 

29,68227 

-38 « 2454 J \ 

11.85444 

-4.10757 

18 

48.99609 

-20.29485  ^ 

xT 

7.17531 

-2.97212 

19 

48.03211 

6.05504  

0.66952 

-3.36596 

20 

26.51649 

26 ,51648  

^FO  - 

;13. 25824 

-32.00824 

21 

-5.97340 

30.03027 

F4 

36.74526 

-11.53291 

22 

-34.64547 

14.35062 

F8 

17.32273 

41.82078 

23 

-46.69306 

-12.7869 7 

F12- 

40.15204 

26.82873 

24 

-0.00000 

37 . 50000  

13,25824 

32.00824 

25 

-6.05503 

48,03210  X. 

^ - 

11.53290 

36.74527 

26 

-20.29482 

48.99606  \ 

17.32272 

41,82075 

27 

-29.68225 

38.24538 

x 

26.82874 

40.15204 

28 

-26.51649 

26.51648^- 

i^-FO 

13.25824 

5.49176 

29 

-17.01077 

25.45843 

F4 

-4.10756 

11.85443 

30 

-14.35062 

34.64544 

F8 

-7.17531 

-2.97210 

31 

-23.97523 

42.05870 

F12 

3.36595 

-0.66952 

FOR  K = 0 » 

8 » 1 6 > 24 

FOR  J = 0 TO  3 

G ( J+K ) » C G ( J+K > +G  < J+K+4 ) 3/2 
G< J+K+4)=F(4*J>*CG<J+K>-G< J+K+4)]/2 


STEP  3 OF  A 32  POINT  FOURIER  TRANSFORM 


FIGURE  5.7 


- 139  - 


* 

OLD  i 

ARRAY 

FACTOR  NEW 

ARRAY 

REAL 

IMAG 

REAL 

IMAG 

0 

12.50000 

0.00000^- 

5^9.37500 

i 0.00000 

1 

10.82532 

0.00000 

5.41266 

0.00000 

2 

6.25000 

0.00000  — = 

~~  ro 

3.12500 

0.00000 

3 

0.00000 

0.00000 

F8 

0.00000 

5.41266 

4 

-0.00000 

0.00000 



«■ — 0.00000 

9.37500 

5 

7.65465 

7.65465 

-3.82733 

11.48198 

6 

0.00000 

18.75000—=: 

— -&^F0 

-0.00000 

-9.37500 

7 

-15.30930 

15.30930 

F8 

3.82733 

11.48198 

8 

18.75000 

-18.75000 

16.00412 

-16.00412 

9 

18.28662 

-15.85990 

ETC. 

12.07262 

-15.00193 

10 

13.25824 

-13.25824 

FO 

2.74588 

-2.74588 

11 

5.85862 

-14.14395 

F8 

0.85798 

6.21400 

12 

18.75000 

18.75000 

2.74588 

2.74588 

13 

-15.85989 

18.28661 

-0.85797 

6.21400 

14 

-13.25824 

-13.25824 

FO 

16.00412 

16.00412 

15 

14.14394 

-5.85861 

F8 

-12.07261 

-15.00192 

16 

13.25824 

-5.49176 

10.21678 

-4.23194 

17 

11.85444 

-4.10757 

6.26198 

-3.73677 

18 

7.17531 

-2.97212 

FO 

3.04147 

-1.25982 

19 

0.66952 

-3.36596 

F8 

0.37080 

5.59246 

20 

-13.25824 

-32.00824 

2.03224 

- 4,90627 

21 

36.74526 

-11.53291 

-1.70339 

7.64791 

22 

17.32273 

41.82078 

FO 

-15,29049 

-36.91451 

23 

-40.15204 

26.82873 

F8 

19. 18082 

38.44865 

24 

-13.25824 

32.00824 

-15.29048 

36.91450 

25 

-11.53290 

36.74527 

-19.18082 

38.44865 

26 

-17.32272 

41.82075 

FO 

2,03224 

-4.90625 

27 

-26.82874 

40.15204 

F8 

1.70339 

7.64792 

28 

13.25824 

5.49176 

3.04147 

1.25983 

29 

-4.10756 

11.85443 

-0.37080 

5.59246 

30 

-7.17531 

-2.97210 

FO 

10,21678 

4.23193 

31 

3.36595 

-0.66952 

F8 

-6.26198 

-3,73675 

FOR  K - 0 

» 4 r 8 » ... 

. « . r 24  r 

28 

FOR  J - 0 r 1 

G ( J+K ) 

= CG ( J+K > + G < J+K+2  > 3/2 

G ( J+K+2 ) -F ( 8* J ) *CG  ( J+K ) -G ( J+K+2 ) 3/2 
STEP  4 OF  A 32  POINT  FOURIER  TRANSFORM 

FIGURE  5.8 


I 


if 


- 140  - 


OLD  ARRAY 


FACTOR 


NEW  ARRAY 


REAL 

IMAG 

REAL 

IMAG 

0 

9.37500 

0.00000 

7 . 39383 

0.00000 

1 

5.41266 

0.00000  "" 

_ 

»-F0 

1.98117 

0.00000 

2 

3.12500 

0.00000  — 

^-1.56250 

2.70633 

3 

0.00000 

5.41266  , 

^►FO 

1.56250 

-2.70633 

4 

-0.00000 

9.37500 

-1.91366 

10.42849 

5 

-3.82733 

11.48198  ETC. 

FO 

1.91366 

-1.05349 

6 

-0.00000 

-9.37500 

1.91366 

1.05349 

7 

3.82733 

11.48198 

FO 

-1.91366 

-10.42849 

8 

16.00412 

-16.00412 

14.03837 

-15.50302 

9 

12.07262 

-15.00193 

FO 

1.96575 

-0.50110 

10 

2.74588 

-2.74588 

1.80193 

1.73406 

11 

0.85798 

6.21400 

FO 

0.94395 

-4.47994 

12 

2.74588 

2.74588 

0.94395 

4.47994 

13 

-0.85797 

6.21400 

FO 

1.80193 

-1.73406 

14 

16.00412 

16.00412 

1.96576 

0.50110 

15 

-12.07261 

-15.00192 

FO 

14.03837 

15.50302 

16 

10.21678 

-4.23194 

8.23938 

-3.98435 

17 

6.26198 

-3.73677 

FO 

1.97740 

-0.24759 

18 

3.04147 

-1.25982 

1.70614 

2.16632 

19 

0.37080 

5.59246 

FO 

.1.33533 

-3.42614 

20 

2.03224 

4.90627 

0.16443 

6.27709 

21 

-1.70339 

7.64791 

FO 

1.86782 

-1.37082 

22 

-15.29049 

-36.91451 

1.94517 

0.76707 

23 

19.18082 

38.44865 

FO 

-17.23566 

-37.68158 

24 

-15.29048 

36.91450 

-17.23565 

37.68157 

25 

-19.18082 

38.44865 

FO 

1.94517 

-0.76708 

26 

2.03224 

-4.90625 

1.86781 

1,37083 

27 

1.70339 

7.64792 

FO 

0.16443 

-6.27709 

28 

3.04147 

1,25983 

1.33533 

3.42614 

29 

-0.37080 

5.59246 

FO 

1.70614 

-2.16631 

30 

10.21678 

4.23193 

1.97740 

0.24759 

31 

-6.26198 

-3.73675 

FO 

8.23938 

3.98434 

FOR  K = 0 t 

2 f 4 t ......  t 

28  » 

30 

_ * 

► 


G(K)  = C:G(K)+G(K+l>3/2 
G(K+1 )=F(0)*CG(K)-G<K+1 )l/2 


# 

OLD  ARRAY 

FACTOR  NEW 

ARRAY 

REAL 

I MAG 

REAL 

IMAG 

0 

7.39383 

0.00000 

7.39383 

0.00000 

1 

1.98117 

0.00000 

8.23938 

-3.98435 

2 

1 . 56250 

2 . 70633s. 

* 14.03837 

-15.50302 

3 

1.56250 

-2.70633  \ 

>/ -17.23565 

37.68157 

4 

-1.91366 

10.42849  v 

S.  >/  -1.91366 

10.42849 

5 

1.91366 

-1.05349 

0.16443 

6.27709 

6 

1.91366 

1.05349  > 

/ \ 0,94395 

4,47994 

7 

-1,91366 

-10,42849 

^ 1.33533 

3.42614 

8 

14.03837 

-15.50302/ 

»1 .56250 

2.70633 

9 

1.96575 

-0.50110 

1.70614 

2.16632 

10 

1.80193 

1.73406 

1.80193 

1.73406 

11 

0.94395 

-4.47994 

1.86781 

1,37083 

12 

0.94395 

4.47994 

1.91366 

1.05349 

13 

1.80193 

-1.73406 

1.94517 

0,76707 

14 

1.96576 

0.50110 

1 , 96576 

0,50110 

15 

14.03837 

15.50302 

1.97740 

0.24759 

16 

8.23938 

-3.98435 

1 .98117 

0,00000 

17 

1.97740 

-0.24759 

1,97740 

-0.24759 

18 

1.70614 

2.16632 

1.96575 

-0.50110 

19 

1.33533 

-3.42614 

J.  .94517 

-0.76708 

20 

0.16443 

6,27709 

1.91366 

-1.05349 

21 

1,86782 

-1.37082 

1.86782 

-1.37082 

22 

1.94517 

0,76707 

1.80193 

-1.73406 

23 

-17.23566 

-37.68158 

1.70614 

-2.16631 

24 

-17.23565 

37.68157 

1.56250 

-2.70633 

25 

1.94517 

-0,76708 

1.33533 

-3.42614 

26 

1.86781 

1.37083 

0,94395 

-4,47094 

27 

0,16443 

-6.27709 

0.16443 

-6 . 27709 

28 

1.33533 

3.42614 

-1.91366 

-10.42849 

29 

1.70614 

-2.16631 

-17.23566 

-37.68158 

30 

1.97740 

0.24759 

14.03837 

15,50302 

31 

8.23938 

3.98434 

8.23938 

3.98434 

BIT  REVERSAL  STEP 
G ( J ) IS  INTERCHANGED  WITH  GOO 
IE  J IS  BIT  REVERSED  K 


I.E.  IF  J * 2 = 00010  BINARY 

32  BIT  REVERSE  J = 01000  = Q DECIMAL 

THEREFORE  G(2)  IS  INTERCHANGED  WITH  GO/ 

STEP  6 OF  A 32  POINT  FOURIER  TRANSFORM 


FIGURE  5.10 


142  - 


CHAPTER  6 : THE  RETRIEVAL  OF  SIGNALS  FROM  NOISE 


6.1  Introduction 

In  many  practical  measurement  situations,  the  signal  under 
observation  is  accompanied  by  significant  amounts  of  "noise"  i.e.,  other 
signals  unrelated  to  the  signal  under  observation  and  adding  no  additional 
information  to  it.  This  noise  may  be  wide  band  (e.g.  Gaussian  or  thermal 
noise),  narrow  band  (e.g.  50  Hz  electrical  mains  interference)  or  a 
mixture  of  both.  The  amplitude  of  the  noise  can  be  sufficient  to  swamp 
the  signal  completely.  In  many  instances,  signal  retrieval  techniques  can 
considerably  improve  the  signal  to  noise  ratio  (SNR) . 

Techniques  such  as  coherence  functions  can  aid  in  the  estimation 
of  the  occurrence  time  of  a signal,  i.e.  the  response  to  a known  stimulus, 
but  where  the  signal  waveform  is  required  other  techniques  must  be  employed. 
This  chapter  contains  a discussion  of  some  signal  retrieval  techniques  and 
their  application..  Further  information  can  be  obtained  from  the 
references  listed  at  the  end  of  the  chapter. 

6.2  The  Retrieval  Process 

It  will  be  evident  from  the  Introduction  that  any  noise  retrieval 
process  when  viewed  in  the  frequency  domain  must  amount,  explicitly  or 
otherwise,  to  a filtering  operation.  Where  the  dominant  components  of  the 
power  spectra  of  signal  and  noise  are  significantly  different,  use  of 
appropriate  high  pass,  low  pass,  band  pass  or  band  stop  filtering  can 
provide  an  adequate  means  of  retrieving  the  signal.  The  realisation  and 
use  of  such  filters  will  be  discussed  later  in  this  Section.  However  in 
many  practical  measurements,  the  signal  and  noise  spectra  overlap  and,  where 
the  signal  is  a response  to  a repetitive  stimulus,  other  techniques  can  prove 
useful;  some  of  these  will  now  be  discussed. 


6.2.1  Coherent  Averaging 

> 

A time-limited  signal  that  is  buried  in  noise  time- locked 
to  an  external  stimulus  can  often  be  retrieved  by  coherent  averaging.  By 
synchronized  additive  combinations  of  repeated  sweeps  of  the  noisy  signal, 
a selective  enhancement  of  the  signal  waveform  can  be  achieved  at  the 
expense  of  the  associated  noise  and  it  can  be  shown  that  this  technique 
produces  a non-linear  filter,  i.e.  one  that  acts  on  the  noise  leaving  the 


signal  unaffected,  even  when  the  two  share  dominant  Fourier  components. 

The  averaging  process  is  achieved  by  synchronising 
each  recorded  sweep  of  data  with  a trigger  waveform  having  a constant 
time  relationship  to  the  signal.  Successive  sweeps  are  averaged,  so 
that,  if  the  signal  repeats  identically  in  form  every  time  it  appears, 
the  averaged  waveform  contains  a constant  amount  of  signal.  If  the 
noise  is  uncorrelated  with  the  signal,  the  summation  of  the  different 
noise  samples  will  reduce  the  noise  by  partial  cancellation  because 
the  contribution  of  the  noise  is  sometimes  positive,  sometimes  negative. 
Thus,  as  will  be  discussed  further  below,  the  noise  is,  in  effect, 
reduced  and  filtered.  The  effect  is  shown  in  Fig.  6.1. 

(a)  Calculation  of  signal  enhancement:  Averaging  earn  be  performed 
using  either  a digital  computer  or  a special-purpose  digital  averager. 

When  writing  a computer  programme,  it  is  convenient  to  store  and  display 
a running  average  of  the  stored  sweeps  of  data.  In  calculating  the  effect 
of  averaging  it  will  be  assumed  that  the  signal  and  noise  are  uncorrelated, 
that  the  noise  is  ergodic  and  that  each  sample  of  it  is  independent  and 
that  the  signal  is  stationary  and  time  invariant. 

If  the  noisy  signal  is  x(t),  then 

x(t)  = s(t)  + n(t)  (6.2.1) 

< 

where  s(t)  is  the  signal  and  n(t)  is  the  noise.  A convenient  measure 
of  the  relative  signal  strength  is  the  " signal- to-noise  ratio",  SNR (t) 
given  by 

SNR(t)  = s(t)  / /(n2)  (6.2.2) 

where  the  denominator  in  eqn  (6.2.2)  is  the  root  mean  square  (RMS)  value 
of  the  noise.  (Note  that  the  signal-to-noise  ratio  is  sometimes  defined 
as  the  ratio  of  signal  and  noise  power.)  It  can  be  shown  that,  after  K 
sweeps  of  data  have  been  accumulated  and  the  results  averaged  the  SNR 
becomes 

I K-l 

SNR(K,t)  = Ks (t)  / /{KR  (o)  + 2 l (K-j)  R (j,t)}  (6.2.3) 

, nn  j=l  ** 


when  R (t)  is  the  auto-correlation  function  of  the  noise.  If  this 
nn 

last  is  known  the  SNR  is  completely  specified.  For  the  case  of 
Gaussian  noise, 


R (0)  = n2 
nn 


R (t)  =0  (t  jt  0) 
nn 


(6.2.4) 


and  (6.2.3)  reduces  to 


SNR  (K, t)  = /K  s(t)  / /(n  ) 


(6.2.5) 


► 


(If  n(t)  is  violently  asymmetrical,  then  the  average  value  will  contain 
a DC  shift  that  has  been  ignored  here) . 

When  n (t)  is  not  Gaussian,  the  averaging  can  lead  to  an 
enhancement  either  much  better  or  much  worse  them  Jv.  depending  on  the 
form  of  the  second  term  in  the  denominator  of  eqn  (6.2.3). 

There  are  two  important  points  to  note  here  : 

(i)  If  n(t)  contains  a significant  oscillatory  term  then 
SNR (K,t)  may  itself  be  an  oscillating  function,  making 
the  use  of  averaging  somewhat  hazardous. 

(ii)  For  a /k  improvement,  it  is  necessary  only  that 


K-l 

£ (K-j)  Rnn<j»t)  = 0 
j=l 


(6.2.6) 


In  addition  to  Gaussian  noise,  this  condition  can  be  shown  to  be  satisfied 
for  any  noise  having  a power  spectral  density  of  the  form  |u>|  ° where 
0 < a (I  (e.g.  "flicker"  noise  in  semi-conductors). 

The  RMS  noise  figures  calculated  above  are  the  mean  or  "expected 
values  of  a stochastic  process.  Thus,  in  practice,  as  the  average  is 
being  built  up,  the  RMS  noise  will  oscillate  about  its  expected  value 
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causing  the  signal  to  rise  out  of  the  noise  somewhat  erratically.  If 

R (t)  is  known,  it  is  possible  to  estimate  the  likely  magnitude  of 
nn 

variability  of  the  noise  fluctuations  in  the  output.  It  is  thus 
possible  to  estimate  confidence  limits  on  any  signal  component  thought 
to  have  arisen  through  the  averaging  process  and  to  separate  stimulus- 
locked  fluctuations  of  the  averaged  noise. 

Where  these  fluctuations  cure  Gaussian  or  near  Gaussian 
applications  of  Student's  t-test  at  the  95%  level  may  be  sufficient. 

Where  the  form  of  the  distribution  is  unknown  confidence  limits  can  be 
estimated  using  the  Chebyshev  inequality  or  by  the  use  of  non-par ame trie 
statistics. 

(b)  Practical  limitations  of  averaging:  In  practice  the  assumptions 
made  in  (a)  above  regarding  the  characteristics  of  the  signal  and  noise 
will  not  all  be  valid.  In  particular, 

(i)  any  latency  "jitter"  in  the  signal  relative  to  its 
stimulus  will  smear  the  averaged  signal, 

(ii)  the  signal  and  noise  will  never  be  entirely  uncorrelated 
since  the  "window"  involved  in  the  sampling  process 
ensures  that  both  have  frequency  components  in  common 
and,  hence,  a finite  cross-correlation.  In  most,  but 
not  all,  instances  this  will  not  be  a significant  factor. 

In  the  absence  of  statistical  checks  the  validity  of  any 
averaged  response  must  always  be  a subject  of  some  concern. 

(c)  Non-linear  filtering:  Since  the  averaging  process  adds  successive 

signal  components  in  phase,  the  frequency  spectrum  of  the  signal  remains 

unaltered.  The  effect  on  the  noise  spectrum,  however,  will  often  be  quite 

considerable.  As  an  example,  consider  the  case  where  sweeps  of  the  signal 

are  taken  at  intervals  of  T seconds.  Then  after  K sweeps  the  total 

o 

recorded  noise  at  any  time  t after  the  commencement  of  a sweep  will  be 


nk  ^ (t)  = n (t)  +H(t  + T ) + n (t  + 2T  )+  •••  +n(t  + (K-l)T  ) 
tot  o o o 


(6.2.7) 


e: 

e 
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On  taking  the  Fourier  transform  of  this  last,  employing  the  shift  theorem, 
and  summing  the  resultant  series  it  will  be  found  that 


N (f ) = N(f)  sin  {KufT  } / sin  {ufT  > 

tot  o o 


(6.2.8) 


i.e.  the  transform  N(f)  has  been  low-pass  filtered.  Thus,  in  this 
instance,  the  averaging  process  acts  as  a non-linear  filter,  filtering 
the  noise  but  not  the  signal. 


(d)  Phase  domain  averaging;  Averaging  in  the  phase  domain  can  be 
used  to  extract  a signal  whose  shape  remains  constant  but  whose  latency 
is  subject  to  considerable  jitter.  The  basis  of  the  technique  is 
illustrated  diagrammatically  in  Fig.  6.2.  It  will  be  seen  that  the  signal 
shape  is  preserved. 


Phase  domain  techniques  are  particularly  effective  in  processing 
noisy  signals  because  most  of  the  information  determining  the  shape  of 
the  signal  is  contained  in  the  phase  spectrum.  The  additions  of  small 
amounts  of  signal  to  noise  will  lead  to  a significant  degree  of  organis- 
ation in  the  phase  spectrum  of  the  noise. 


(e)  A word  of  warning:  In  the  practical  realisation  of  data  acquisition 
and  averaging,  due  consideration  should  be  given  to  the  choice  of  word 
length.  If,  for  example,  the  signal  of  interest  is  immersed  in  40  db  of 
noise,  the  6 most  significant  bits  of  the  word  will  contain  noise  only 
and  it  will  be  seen  that  the  use  of  10-bit  A-to-D  conversion  will 
considerably  degrade  the  signal- to-noise  ratio  of  the  acquired  signal  and 
prolong  the  averaging  process. 


The  effects  of  rounding  and  truncation  errors  on  the  addition  and 
division  operations  needed  to  compute  the  averaged  signal  can  be  quite 
significant  if  double  precision  arithmetic  is  not  used. 


6.3  Digital  Filtering 


6.3.1  General 


In  many  instances,  averaging  or  other  non-linear  signal 
processing  techniques  are  either  inapplicable  or  yield  unsatisfactory 
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results.  Where  the  residual  signal  and  noise  have  significantly 
different  frequency  spectra,  linear  filtering  may  effect  a significant 
improvement  in  the  signal-to-noise  ratio.  Although  the  specification 
of  a suitable  filter  clearly  implies  some  "a  priori"  knowledge  of 
signal  and  noise  spectra,  in  most  cases  enough  information  is  avail- 
able to  allow  a reasonable  specification  to  be  drawn  up  initially  even 
if  later  changes  have  to  be  made. 

Where  the  signals  are  already  stored  in  a digital 
computer  the  use  of  digital  techniques  enables  extremely  precise  filter 
characteristics  to  be  realised,  as  well  as  saving  the  costs  involved  in 
the  construction  of  complex  analogue  filters.  In  addition  digital 
filters  can  provide  characteristics  unrealisable  in  analogue  form  (with- 
out the  use  of  delay  times) . 

Digital  filtering  can  be  produced  by  time  or  frequency 
domain  manipulations  of  data  and  these  two  approaches  will  now  be  consid- 
ered in  detail. 

6.3.2  Time  Domain  Filters 

Time  domain  averages  are  realised  using  moving  average 
techniques  analogous  to  the  moving  average  series  discussed  in  Chapter  4. 
The  sampling  process  converts  the  signal  time  series  {x(t)}  to  the 
stored  data  array  {x(n)}  by  sanpling  it  at  intervals  of  At  sec  (where 
the  interval  is  sufficiently  small  to  avoid  aliasing  problems) . 


(a)  Non- recursive  filters:  A non-recursive  filter  merely  replaces 

x(n)  with  a weighted  sum  y(n)  made  up  from  x(n)  and  neighbouring 
points  in  the  data  array.  This  may  be  expressed  in  terms  of  a difference 
equation 


a^  x(n-j) 


(6.3.1) 


where  m,  r * 0,  1,  2 ... 

Such  a filter's  operation  may  be  illustrated  by  the 
symmetrical  three  point  filter  series 


3mm 


(Note  that  this  filter  will  be  physically  unrealisable  since  x(n+l) 
is  ahead  of  x(n)  in  time.  This  is  of  no  concern  where  X is 
already  stored.) 


Since  the  sampling  rate  is  l/(At)  the  filter  transfer 


function 


(where  w = 2irf  and  X and  Y are  the  Fourier  transforms  of  x and  y) 
may  be  found  by  noting  that,  from  the  shift  theorem , the  transform  of 
x(n-j)  is 


X (oj)  exp(-ijo»(At) ) 


Y(u>)  = (1/2)  X(a>)  + (1/4) {exp (ito (At) ) + exp (-ico (At) ) }x (w) 


This  transfer  function  is  shown  in  Fig.  6.3.  Note  that  the  data  band 
width  is  limited  by  the  sampling  rate  to 


ir/2(At)  rad/sec 


so  that  only  the  first  "quadrant"  of  the  filter  expression  is  of 
significance.  Note  also  that  the  low  pass  filter  introduces  zero  phase 
shift.  This  is  a result  of  the  symmetrical  filter  algorithm  and  is  a 
useful  property  where  measurement  of  the  latency  of  components  in  the 
filtered  signal  is  important.  The  effect  of  a filter  containing  a pair 
of  critically  damped  conjugate  poles  is  included  for  comparison. 
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The  corresponding  high-pass  filter 
H (w)  * (1  - cos  ojt)  /2  (6.3.7) 

is  obtained  using  the  series 

y (n)  = -(1/4)  x(n-l)  + (1/2)  x(n)  - (1/4)  x(n+l)  (6.3.8) 

These  filters  can  be  shown  to  have  3 db  cut-off  points  of 

f = 0.183  / (At)  Hz  (6.3.9) 

o 

The  filtering  operation  may  be  repeated  as  many  times  as  desired. 

This  has  the  effect  of  lowering  f in  the  low-pass  case,  or  raising 
f in  the  high  pass  case  by  a factor 

(l-l/2k) 

0.875  arcos  {2  - 1}  (6.3.10) 

after  k operations. 

The  design  of  digital  filters  having  given  transfer 
functions  can  be  carried  out  in  the  z plane  using  the  conformal 
transformation 

z = exp  (iut)  (6.3.11) 

This  permits  the  reduction  of  the  difference  euqation  to  a transfer 
function  e.g.  for  the  above  filter  (6.3.2)  use  of  (6.3.11)  gives 

Y (z)  = (1/4)  (1/z  + 2 + z)  X (z)  (6.3.12) 

i.e.,  H (z)  = {1  + (z  + 1/z)  /2}  /2  (6.3.13) 

Discussion  of  z plane  design  techniques  can  be  found  in  various  texts. 


(b)  Recursive  (auto-regressive)  filters:  While  the  non-recursive 
filter  is  extremely  simple  the  algorithms  become  large  and  unwieldy  where 
high  performance  filters  are  required.  In  these  situations,  recursive 
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algorithms  (analogous  to  the  autoregressive  series  of  Chapter  4)  may 
be  of  use.  These  include  previously  filtered  data  points  in  each 
moving  average  calculation.  This  leads  to  a drastic  reduction  in 
the  number  of  multiplications  required. 

The  general  difference  equation  may  be  written  as 


y (n)  = E a.x(n-j)  + E b y(n-j) 
j=-m  ^ j=l-m 


(6.3.14) 


The  power  of  recursive  filters  is  best  illustrated  by  a simple  example. 
Consider  the  non-recursive  series 


y (n)  = E x(n+j) 
j=-P 


(6.3.15) 


comprising  (2p+l)  unit  samples,  or  transforming  into  the  z plane 


H(z)  = Y (z)  / X(z)  = (z-p  + •••  + 1 + •••  + z?) 
L.e. , H (z)  = {z5  - z (p+1)}  / {1  - z 1} 


On  applying  the  inverse  transformation,  the  result  is 


y (n)  - y(n-l)  = x(n+p)  - x(n-(p+l)) 


i.e.,  y (n)  = y(n-l)  + x(n+p)  -x(n-(p+l) 


(6.3.16) 


(6.3.17) 


(6.3.18) 


(6.3.19) 


This  recursive  difference  equation  contains  three  points  only.  Note 
that  this  filter  algorithm  is  equivalent  to  a non-recursive  algorithm 
containing  any  number  of  points. 

Recursive  filters  may  need  some  additional  non-recursive  filters 
to  remove  considerable  side  lobes.  Coefficient  rounding  and  truncation 
can  lead  to  instability  and  oscillation.  With  these  qualifications 
recursive  filters  can  provide  a useful  tool  in  signal  processing.  They 
are  particularly  attractive  for  the  design  of  real  time  digital  filters 
utilising  shift  registers  and  digital  multipliers. 
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Although  only  high  and  low  pass  filters  have  been 
discussed  above,  it  will  be  evident  that  band-pass  and  band-stop 
filters  can  be  realised  using  the  same  techniques. 


6.3.3  Fourier  Transform  Filters 


Transformed  data  may  be  manipulated  in  the  frequency 
domain  by  multiplying  the  real  and  imaginary  components  by  the 
appropriate  weighting  function.  The  modified  data  may  then  be 
returned  to  the  time  domain.  Very  high  performance  filters  can  be 
produced  using  this  technique  which  is,  in  theory,  capable  of 
producing  "brick  wall  filters"  having  an  infinite  asymptotic  slope  at 
the  cut-off  frequencies. 

In  practice,  the  complete  removal  of  one  or  two  spectral 
components  as  a means  of  removing  a single  interfering  signal  can  lead 
to  unacceptable  "ringing"  in  the  filtered  signal.  This  effect  can  be 
reduced  by  tapering  the  spectral  conponents  each  side  of  the  components 
being  removed.  Where  only  simple  filters  are  required  time  domain 
filtering  will  in  general  be  simpler  and  more  convenient. 

4 Selection  of  Filter  Characteristics  - Optimal  Filterin 


Choice  of  a filter  to  separate  signal  and  noise  is 
normally  approached  on  a somewhat  Empirical  basis.  Where  there  is  a 
priori  knowledge  of  the  spectral  characteristics  of  signal  and  noise  it 
is  possible  to  arrive  at  an  optimal  filter  characteristic  for  signal 


recovery. 


Filter  characteristics  generated  by  this  technique  are 
often  physically  unrealisable  but  this  is  unimportant  where  digital  filters 
are  used.  The  filters  discussed  here  will  be  assumed  time-invariant,  but 
the  same  approach  may  be  extended  to  signal-dependent  (e.g.  Kalman) filters  • 
Consider  a signal  s(t),  having  a power  spectral  density  Pg(f),  immersed 
in  noise  n(t),  having  a power  spectral  density  Pn<f).  After  signal 
and  noise  are  passed  through  a filter  H(f)  (the  inverse  of  which  is  h(t)), 
the  signal-to-noise  ratio  is 


't 
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SNR(t)  = {/^  s(t-X)  h^(X)  dX}  / {/^n  (t-X)  h (X)  dX}  (6.3.20) 


Choosing  H(f)  to  maximise  SNR  leads  to  the  result 
H (f ) = S(f) 


(6.3.21) 


i.e.,  the  optimal  filter  should  follow  the  a priori  estimate  of  the 
signal  spectrum  S (the  transform  of  s ) . 

Another  approach  is  to  find  the  filter  that  makes  the 
output  s + n,  the  best  least  squares  estimate  of  the  signal.  This 
leads  to  the  result 

H(f)  = P (f)  / P (f)  (6.3.22) 

s n 

This  operation  will  practically  always  lead  to  a physically  unrealis- 
able  filter.  In  most  cases  the  simple  result  (6.3.21)  will  be  adequate. 


6.4  Spectral  Analysis 

6.4.1  General 

In  many  instances  the  signal  being  measured  is  continuous 
(e.g.  mechanical  vibration  in  response  to  an  applied  force).  In  such 
cases  the  signal  is  best  characterised  by  its  frequency  spectrum.  This 
normally  takes  the  form  of  a power  spectrum,  indicating  the  distribution 
of  energy  with  frequency.  (This  has  the  advantage  of  being  easier  to 
interpret  and  use  than  the  magnitude  spectrum  in  many  engineering 
situations.) 


Any  finite  data  sample  can  supply  only  an  estimate  of  the 
true  power  spectrum?  thus  spectral  computation  is  meaningless  without  an 
assessment  of  the  relevant  confidence  limits. 

6.4.2  Power  Spectrum  Estimation 

Once  the  spectral  resolution  and  confidence  limits  required, 
and  the  number  of  data  points  in  each  FFT  are  known,  it  is  possible  to 
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specify  the  experimental  protocol.  This  has  already  been  discussed 
in  Chapter  5 but,  for  completeness , will  be  outlined  again  here. 


The  sampling  interval,  At,  is  of  course  governed  by 
the  Nyquist  criterion 


At  4 1/2  f 


max 


(6.4.1) 


where  fffny  is  the  highest  frequency  component  in  the  signal.  The 


number  of  points,  N,  in  each  data  block  is  linked  to  the  frequency 
resolution,  Af  by 


Af  = 1/ (N  At) 


(6.4.2) 


It  thus  pays  to  use  as  large  a value  of  N as  possible.  The  total 
length  of  signal  required  is  determined  by  the  amplitude  uncertainty 
that  can  be  tolerated  in  the  computed  spectrum.  Provided  that  the 
spectrum  is  reasonably  flat  or  has  been  "pre-whitened"  by  the  use  of 
a suitable  filter  then  adjacent  spectral  harmonics  will  be  essentially 
uncorrelated  and  will  follow  a chi-squared  distribution  having  2 degrees 
of  freedom  per  harmonic  (since  each  contains  a real  and  an  imaginary 
contribution) . The  chi-squared  distribution  having  k degrees  of 
freedom  has  a mean  k and  a variance  2k.  Hence,  the  coefficient  of 
variation,  CV,  is 


CV  = standard  deviation  / mean  = V(  2/k) 


(6.4.3) 


In  particular,  if  k = 2,  CV  = 1,  and,  however  long  the  data  run,  the 
computed  spectra  will  be  unreliable.  The  problem  cam  be  overcome  either 
by  the  use  of  frequency  smoothing  or  segment  averaging  to  increase  the 
effective  number  of  degrees  of  freedom. 


(a)  Frequency  smoothing;  If  m neighbouring  values  of  the  computed 
spectrum  are  averaged  then  the  resolution  Af  will  be  degraded  to 
m(Af)  but  each  point  of  the  new  spectrum  will  have  2m  degrees  of 
freedom.  The  spectrum  will  therefore  have  a CV  given  by 


CV  = l//m 


(6.4.3) 
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Note  that,  using  this  method,  a total  data  length 

Tr  = mN (At)  (6.4.4) 

will  be  required  to  maintain  the  original  resolution. 

(b)  Segment  averaging:  An  alternative  approach  is  to  average 
corresponding  values  from  q separate  spectra.  This  will  lead  to  a 
CV  given  by 

CV  = l//q  (6.4.5) 

and  a total  data  run  of 

= qN (At)  (6.4.6) 

The  above  two  techniques  may  be  combined  if  desired.  Thus  for  a given 
value  of  CV,  it  is  possible  to  specify  completely  a procedure  for 
the  calculation  of  power  spectra.  Tables  of  the  chi-squared  distribut- 
ion for  k degrees  of  freedom  are  readily  available,  but  many  of  these 
do  not  include  the  large  values  of  k often  required  for  data  processing. 
For  k > 50,  the  relationship 

X2  — { z + /(2k  - 1)}  /2  (6.4.7) 

P P 

where  z is  the  corresponding  percentile  deviation  for  a normal 
P 

distribution  will  often  be  found  useful. 

(c)  Effective  degrees  of  freedom:  The  final  spectrum  can  be  used  to 
calculate  the  number  of  effective  degrees  of  freedom,  kg.  This  last  will 
be  less  than  k since,  in  general,  the  assumption  of  independence  of 
adjacent  spectral  values  will  not  be  entirely  valid.  The  equivalent  band- 
width B of  the  power  spectrum  P(f)  is  given  by 

E 

b = (7p(f)  df}2  / ;"{p(f)}2  df 
E o o 


(6.4.8) 


In  terms  of  this 


k = 2 B„  T - 1 (6.4.9) 

e E r 

This  result  may  be  used  to  correct  the  confidence  limits  of  the 
computed  spectrum. 

(d)  Bias  in  power  spectral  estimation:  Any  estimate  of  a process 
based  on  a finite  amount  of  information  will  generally  display  some 
bias,  however  small.  This  may  be  defined  by 

b = G - G (6.4.10) 

e 

where  b is  the  bias 

G is  the  true  value  of  the  function 

G is  its  estimated  value, 
e 

In  ref.l  it  is  shown  that  the  normalised  mean  square  error  of  a power 
spectrum  is 

t2  = (l/bTr)  + (b4/576)  (PU/P)2  (6.4.11) 

The  second  term  represents  a bias  which  becomes  important  for  a spectrum 
having  significant  dips  or  peaks.  It  will  thus  be  evident  that  an 
extremely  careful  choice  of  b must  be  made  when  computing  such  spectra. 
It  earn  be  shown  that  the  effect  of  this  bias  error  is  to  compress  the 
dynamic  range  of  the  spectrum,  minimising  the  amplitudes  of  peaks  or 
troughs . 
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7.1  Introduction 

Fourier  transforms  are  a powerful  tool  for  solving  many  vibration 
problems  and  their  use  is  finding  more  and  more  acceptance  each  day. 

Their  use  allows  a more  rapid  determination  of  the  natural  modes  and 
frequencies  of  a structure  without  a consequent  loss  of  accuracy.  This 
forms  the  basis  of  the  so-called  "chirp"  test  wherein  the  structure  is 
excited  by  rapid  frequency  sweeps  and  the  transient  response  measured  at 
various  locations  on  the  structure.  The  required  natural  modes  and 
frequencies  are  then  calculated  from  the  measured  data,  contrasting 
with  the  so-called  "steady- state"  test  in  which  each  of  the  modes  is 
excited  in  turn  and  individually  measured. 

The  most  important  concept  in  this  approach  to  vibration  analysis 
is  that  of  a transfer  function  of  a system.  This  function  is  a mathematical 
description  of  the  system  and  describes  the  way  in  which  it  responds  to  an 
input.  If  the  system  is  a mechanical  one,  i.e.  if  it  consists  of  masses, 
springs  and  dampers  then  the  transfer  function  describes  the  relationship 
between  an  exciting  force  and  the  resulting  vibration.  Similarly  the 
transfer  function  of  an  electrical  system  may  relate  a current  at  some 
point  in  the  circuit  to  an  applied  voltage. 

To  help  illustrate  the  properties  of  a transfer  function  the 
equations  of  motion  of  a simple,  single-degree-of-freedom,  mechanical 
system  will  be  developed  in  Section  7.2.  The  extension  to  a multi-degree- 
of-freedom  system  will  be  indicated  but  not  rigorously  performed.  The 
following  Sections  then  describe  some  of  the  experimental  methods  used  to 
determine  the  transfer  function  . 

7.2  Linear  System  Theory 

Consider  the  simple  mass-spring-damper  system  shown  in  Fig. 7.1. 

If  the  mass  is  excited  by  a force  f (t) , it  responds  such  that  its 
amplitude  of  vibration  is  given  by  x(t).  This  simple  system  possesses 
two  important  properties  which  are  fundamental  to  the  application  of 
Fourier  transform  theory.  Firstly  it  is  linear.  That  is,  if  a force 
f (t)  produces  a response  x (t) , and  another  force  f_(t)  produces  a 
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different  response  x2(t) , then  a combined  force  f 1 (t)  + f2(t)  will 

produce  a response  x^(t)  + x2(t).  This  implies  that  a complicated  # 

forcing  function  may  be  considered  in  terms  of  simpler  constituent 

functions,  the  response  to  each  determined,  and  then  these  responses 

combined  to  give  the  response  to  the  original  force.  The  second  « 

important  property  is  that  of  shift  invariance.  This  implies  that  if 

a force  is  shifted  in  time  (i.e.,  delayed  or  advanced)  the  only  effect 

is  that  the  response  is  shifted  by  the  same  amount. 

The  second-order,  linear  differential  equation  relating  the 
response  x(t)  to  the  applied  force  f(t)  is 

m x (t)  + b x (t)  + k x (t)  = f(t)  (7.2.1) 

where  m is  the  mass 

b is  the  viscous  damping  coefficient  1 

and  k is  the  stiffness. 

« 

Given  a particular  f(t),  this  eqn  (7.2.1)  may  then  be  solved  for  x(t) 
by  the  standard  mathematical  techniques. 

For  our  purposes  here  we  will  consider  the  special  case  of  a simple 
harmonic  forcing  function.  This  may  be  written  a number  of  ways  but  it  is 
easiest  to  use  the  complex  exponential  notation  : 

f(t)  = F.exp(iwt)  = F. (coswt  + i sin  wt)  (7.2.2) 

where,  as  is  usual,  only  the  real  or  the  imaginary  part  represents  the 
relevant  physical  quantity.  Here  F is  the  peak  amplitude  of  the  force 
which  oscillates  at  a circular  frequency  of  a>  radians/second.  Since  the 
system  is  linear  the  steady  (non-transient)  response  x(t)  will  also  be 
simple  harmonic  of  frequency  <u  but  with  a phase  lag  of  $.  This  response 
function  may  also  be  written  in  a number  of  ways,  such  as 

* 

x (t)  - A exp(ioit  + i$)  - (V  + i Z)  exp  (iwt)  - X exp(iwt) 

(7.2.3) 
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where  A,  Y and  Z are  real  but  X is  complex . He  will  use  the  last 


form,  i.e. 


x(t)  = X exp  (ib>t) 


Substituting  eqns  (7.2.2)  and  (7.2.4)  into  eqn  (7.2.1)  gives 


(-to  m + itob  + k)  X exp  (itot)  - F.exp  (itot) 


i.e.. 


(7.2.4) 


(7.2.5) 


X = F / (-to  m + iiob  + k) 
* F H (to) 


(7.2.6) 


H(u)  is,  of  course,  complex  and  its  modulus  and  argument  (or  phase)  are 
given  by 


|h|  = l//((-to2  m + k)2  + to2  b2) 


arg(K)  = arctan  ((-tob)/(-to  m + k) ) 


(7.2.7) 


(7.2.8) 


The  relationships  given  in  eqns  (7.2.7)  and  (7.2.8)  are  illustrated  as 
functions  of  to  in  Figs.  7.2a  and  7.2b.  The  two  figures  may  be  combined 
to  give  a single  vector  plot  or  Argand  diagram  as  shown  in  Fig.  7.2c. 

It  may  be  worth  noting  here  that  since  F is  real  the  argument 
of  H is  the  same  phase  lag  4>  as  given  in  eqn  (7.2.3).  The  natural 
frequency  of  the  system  is  defined  to  be  that  frequency  for  which  this 
phase  is  -II/2  radians  and  Fig.  7.2b  shows  that  this  frequency  is 
/ (k/m)  radians/second.  The  peak  of  the  modulus  curve  shown  in  Fig.  7.2a 
occurs  at  a lower  frequency,  / ((k/m)  - (b2/2m2)).  The  height  of  the  peak 
in  Fig. 7. 2a  is  given  by 


| H I * 2m/(b/(4mk  - b2)) 

n 


(7.2.9) 


and  the  so-called  "static  deflection"  (i.e.,  |h|  at  aero  frequency)  is 


I H (o)  I = 1A 


(7.2.10) 


9 ■■■■■I 
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Then  the  magnification  at  resonance  is  given  by 


R = |h|m  / |H(o)|  = 2m  k / ( b /(4m  k - b )) 


(7.2.11) 


The  resonant  frequency  (or  peak  amplitude  frequency)  is  given  by 


“m  = “ fa2/(2m2)) 


If  the  values  of  R and  u cure  known  the  two  parameters 

M 


(7.2.12) 


u = / (k/m)  = natural  frequency 


(7.2.13) 


and  y = b/(2  /(mk))  = ratio  of  damping  to  critical 


(7.2.14) 


may  be  evaluated  exactly. 


If  the  damping  is  small  (y  { 0.05,  i.e.  4 5%  of  critical  damping)  it  can 
be  shown  that  the  following  approximate  relations  may  be  used  : 


w 5 (1  + 1 / (4  R ))  u> 


(7.2.15) 


and  y s 1/(2R) 


(7.2.16) 


These  two  important  parameters  which  characterise  the  system  may  also  be 
obtained  from  the  vector  plot  illustrated  in  Fig.”7. 2c.  Again  if  the 
damping  is  small  (y  4 0.05)  this  plot  is  almost  circular  with  the  negative 
imaginary  axis  as  a diameter.  The  natural  frequency,  u,  is  the 
frequency  at  which  this  diameter  intersects  the  circle.  Then  if  a second 
diameter  drawn  normal  to  the  first  intersects  the  circle  at  frequencies 

u>_  and  ai 

A B 


Y = (u>B  - uA)  / (2  u>  ) 


(7.2.17) 


Whilst  this  derivation  of  the  form  of  H(o))  is  for  a single  degree-of-freedom 
system,  the  extension  to  many  degrees-of-freedom  is  straight-forward. 
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Eqn  (7.2.5)  then  becomes  a set  of  simultaneous  equations,  usually 
written  in  matrix  torm.  The  exact  form  of  H(w)  is  then  more  difficult 
to  write  down  but  the  plot  of  its  modulus  will  be  similar  to  that  for 
the  single  degree-of-freedom  system  (Fig.  7.2a)  except  that  there  will, 
in  general,  be  more  peaks  (Fig.  7.3).  Each  peak  indicates  a resonance, 
and  there  is  a natural  frequency  u>,  and  damping  ratio  y,  associated 
with  each  resonance.  If  the  peaks  are  well  separated,  each  may  be 
considered  in  isolation,  and  the  single-degree-of- freedom  eqns  (7.2.15) 
and  (7.2.16)  may  be  applied.  If  however,  some  of  the  peaks  are  grouped 
closely  together,  more  sophisticated  techniques  must  be  used. 

7.3  Transfer  Function 

The  function  H(u)  which  relates  the  simple  harmonic  response  X, 
to  the  simple  harmonic  applied  force  F,  is  the  transfer  function  of  the 
system.  This  transfer  function  is  rarely  known  explicitly  but  must  be 
determined  experimentally,  usually  as  a table  of  complex  values  at 
discrete  frequencies  over  some  range  of  interest.  It  may  be  determined 
directly  by  applying  a simple  harmonic  force  of  some  frequency  u to  the 
system  and,  after  all  the  transients  have  decayed,  measuring  the  simple 
harmonic  response.  This  gives  the  value  of  H at  the  one  frequency  u>. 

The  frequency  is  then  altered  and  the  process  repeated  until  the  whole 
of  Fig.  7.3  is  built  up.  This  is  a very  slow  process  and  it  is  far 
quicker,  theoretically  at  least,  to  use  Fourier  transforms. 

Consider  any  arbitrary  forcing  function  f (t) . This  may  be  expressed 
in  terms  of  its  Fourier  integral  as 

f(t)  = (1/2*)  /”f(u>)  exp  (iwt)  du>  (7.3.1) 

This  implies  that  the  arbitrary  function  f (t)  may  be  broken  down  into  an 
integral  sum  of  simple  harmonic  components  over  a continuous  range  of 
frequencies.  Then  each  component  produces  a response  given  by  eqn  (7.2.6) 
as 

X(u>)  = H (b>)  F(oi)  ' (7.3.2) 

Due  to  the  linearity  and  shift- invariance  of  our  system  the  total  response 
x(t)  is  the  integrated  sum  of  the  component  responses,  viz 


x(t)  * (l/2x)  /"  X(w)  exp  (iwt)  dw  (7.3.3) 

wOO 

Then  from  eqn  (7.3.2)  the  transfer  function  is  given  by 

H(<d)  = X(w)  / F (w)  (7.3.4) 

where,  as  shown  in  the  first  lecture,  the  frequency  domain  functions  are 
given  by 

X(co)  = F x(t)  exp (-iwt)  dt  (7.3.5) 

•09 

and 

F(w)  = F f(t)  exp  (-iwt)  dt  (7.3.6) 

•pQO 

The  eqn  (7.3.4)  may  be  considered  as  the  definition  of  the  transfer 
function  H (w) . 

Referring  again  to  eqn  (7.3.2)  it  can  be  seen  that  this  is  in  the 
form  of  a product  and  so,  from  the  convolution  theorem, 

x(t)  = h(t)  * f(t)  (7.3.7) 

where  * indicates  convolution  and  where  h(t)  is  the  Fourier  transform 
of  H(w) . Consider  the  special  case  of  f(t)  = 6(t)  where  6 (t)  is  the 
delta  function,  i.e.  an  impulse.  Then  the  associated  response  x^t) 
is  given  by 

x (t)  = h(t)  * 6 (t) 

= h(t)  (7.3.8) 

Thus  h(t) , the  Fourier  transform  of  the  transfer  function,  is  the  impulse 
response  of  the  system  and  is  given  by 

h(t)  * (1/2it)  /°°H(u>)  exp  (iwt)  dw  (7.3.9) 

—00 

For  the  single  degree-of- freedom  system  described  earlier  this  gives 


h(t)  = A exp  (-At)  sin  (w  t) 


(7.3.10) 


where  A is  an  amplitude  factor  = Wp/m 

wD  is  the  circular  decay  frequency  * to  /(I  - y ) and  X is  the 
decay  rate  » yu).  A typical  plot  of  h(t)  is  shown  in  Fig.  7.4.  For  a 
multi-degree-of- freedom  system  the  form  of  the  impulse  response 
involves  combinations  of  terms  like  the  one  in  eqn  (7.3.10). 

The  physical  interpretation  of  eqn  (7.3.2)  is  straightforward. 

The  forcing  function  f(t)  is  considered  to  be  the  sum  of  a number  of 
harmonic  oscillations,  the  amplitudes  of  which  cure  specified  by  F(o>). 

The  transfer  function  H(w)  indicates  how  the  system  responds  to  each 
frequency,  giving  the  responses  as  X (w) . The  equivalent  equation  in  the 
time  domain  is  eqn  (7.3.7).  The  physical  significance  of  this  equation 
is  as  follows.  The  forcing  function  f (t)  is  considered  to  be  a train  of 
impulses  of  the  appropriate  amplitude  and  time  delay.  The  response  of 
the  system  to  each  of  these  impulses  is  obtained  and  the  convolution 
integral  sums  them  to  give  the  final  response. 

(Whilst  in  the  above  the  Fourier  transform  has  always  been  written 
in  the  integral  form,  in  practical  calculations  it  generally  proves 
necessary  to  approximate  this  by  a DFT  which,  of  course,  is  evaluated 
using  the  FFT  algorithm.  This  should  be  borne  in  mind  throughout.) 

7.4  Determination  of  the  Transfer  Function 

It  was  shown  in  eqn  (7.3.4)  that  the  transfer  function  could  be 
determined  from  the  Fourier  transforms  of  the  input  force  and  of  the 
resulting  response  of  the  system.  This  is  true  for  any  forcing  function 
and  we  will  first  consider  the  special  case  of  the  impulsive  force. 

Eqn  (7.3.9)  shows  that  the  impulse  response  of  the  system  is  the 
Fourier  transform  of  the  transfer  function.  This  then  indicates  an 
apparently  easy  method  of  determining  the  transfer  function  - i.e.,  by 
applying  an  impulse  and  then  taking  the  Fourier  transform  of  the  response. 
The  difficulty  in  this  approach  is  in  generating  the  impulsive  force.  It 
must  be  of  infinite  amplitude  and  zero  duration,  a pair  of  conditions 
impossible  to  meet.  Even  an  approximation  to  this,  e.g.  a hammer  blow 
for  mechanical  systems  or  a voltage  pulse  for  electrical  systems,  is  not 
always  successful  as  it  requires  a high  energy  concentration  (i.e.  large 


amplitude)  in  order  to  generate  sufficient  length  of  time-history 
response  for  accurate  analysis*  especially  if  the  damping  (b)  is  high. 

This  application  of  a high  energy-concentration  pulse  is  not  always 
acceptable  as  it  may  damage  the  system  under  investigation.  It  is 
better  if  the  input  energy  is  applied  at  a lower  level  to  the  system 
but  for  a longer  time. 

In  this  latter  approach  the  determination  of  the  transfer 
function  is  achieved  using  eqn  (7.3.4)  directly  with  a convenient  form 
of  input  force.  The  calculation  is  more  accurate  if  the  applied  force 
has  nearly  equal  energy  at  all  frequencies  over  the  range  of  interest 
and  there  are  three  common  forcing  functions  with  this  property.  The 
first  is  white  noise  which  is  defined  as  a random  signal  which  has 
equal  energy  at  all  frequencies,  provided  that  the  noise  is  of  sufficiently 
long  duration.  Usually  it  cannot  be  obtained  in  the  real  world  so  a pseudo- 
white noise  is  used.  This  signal  has  nearly  equal  energies  over  some 
limited  range  of  interest.  This  forcing  function  and  its  power  spectrum  are 
shown  in  Fig.  7.5a. 

A more  common  forcing  function  is  the  fast  frequency  sweep.  If  the 
transfer  function  is  required  over  scane  range  to  u>^,  then  the  frequency 

of  the  force  is  increased  rapidly  from  to  whilst  exciting  the 

system.  This  force  and  its  power  spectrum  are  shown  in  Fig.  7.5b.  This  is 
perhaps  the  most  comnonly  used  force  input. 

The  third  forcing  function  is  called  the  "impulse  sine"  and  is  the 
same  as  the  sine  function  described  in  the  first  lecture.  It  is  of  the  form 

f(t)  = sin (u1 (t-T) ) / (t-T)  (7.4.1) 

This  function  and  its  power  spectrum  are  shown  in  Fig.  7.5c.  The  power 
spectrum  is  essentially  flat  up  to  the  cut-off  frequency  and  in  the 
limit  as  approaches  infinity,  the  impulse  sine  becomes  the  delta  function 
However,  it  is  easier  to  generate  than  a true  impulse  and  it  has  a lower 
energy  concentration.  Moreover  it  applies  energy  to  the  system  over  a 
relatively  wide  frequency  range  much  more  quickly  than  a' frequency  sweep. 

For  these  reasons  it  is  a very  useful  forcing  function  especially  when 
test  conditions  cannot  be  maintained  for  more  than  a brief  period. 
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The  foregoing  method  of  determining  the  transfer  function  involved 
dividing  the  Fourier  transform  of  the  response  by  that  of  the  input 
force.  However  the  transfer  function  may  also  be  calculated  from  a 
knowledge  of  the  cross-  , and  auto-power  spectra  of  the  response  and 
input. 


The  auto-power  spectrum  of  the  response  X (w)  was  defined  in 
earlier  lectures  as 

GXX(W)  = X(U>)  X*  (a>)  (7.4.2) 

where  * indicates  the  complex  conjugate.  Then  substituting  from  eqn 
(7.3.2)  gives 


GXXfa>)  = (H(“)-F  <“>>)•  («(“)•  F (“>))* 

= |h(u>)  |2.  Gtc (co)  (7.4.3) 

FF 

where  G (u)  is  the  auto-power  spectrum  of  the  exciting  force.  Thus 
FF 

|h(u>)  I = G^(w)  /Gpp (m)  (7.4.4) 

Eqn  (7.4.4)  allows  the  calculation  of  only  the  modulus  of  the  transfer 
function  with  no  phase  information.  However  this  is  often  sufficient  to 
determine  many  of  the  properties  (e.g.,  natural  frequencies  and  dampings) 
of  the  system.  Eqn  (7.4.4)  may  also  be  used  when  the  exact  form  of  the 
input  is  unknown  but  where  its  power  spectrum  is  known  as  would  be  the 
case  if  the  excitation  were  electrical  noise.  For  this  case  the  energy 
spectrum  of  the  input  is  frequently  assumed  to  be  flat  over  the  frequency 
range  of  interest. 

In  order  to  obtain  the  phase  content  of  the  transfer  function  from 
the  power  spectra,  the  cross-power  spectrum  must  be  used.  This  is  defined 
as 

G (u>)  = X(«).  F*  (w) 

XF  » 

= (H(u).  F (u>) ) . F*  (w) 


H (<*>)•  (o>) 


(7.4.5) 


whence 


H (w)  = Gyglu)  / G^to)  (7.4.6) 

Note  that  eqn  (7.4.6)  may  only  be  applied  if  the  forcing  function  f(t) 
is  known.  The  cross-power  spectrum  may  then  be  determined  from  the 
cross-correlation  function  or  from  the  Fourier  transforms  of  the  input 
and  response. 

7.5  Effect  of  Noise 

In  the  preceding  Sections  it  was  assumed  that  the  Fourier  transform 
of  both  the  input  force  and  of  the  response  were  obtainable  and  also  that 
they  were  both  precise  and  accurate.  If,  as  usually  occurs  in  the  real 
world,  the  signals  are  contaminated  with  noise  then  ways  to  overcome  this 
must  be  found. 

Consider  now  the  case  where  the  response  X(u>)  is  partially  due  to 
the  input  force  F(u>),  and  partially  due  to  some  random  noise  N(w). 

i.e. , 

X(a>)  = H(u>) . F (to)  + N(u>)  (7.5.1) 

Then  if  eqn  (7.3.4)  is  used  to  determine  the  noise  contaminated  transfer 
function  (u>) , we  get 

H^io)  = X(u>)  / F (to) 

= H(u>)  + N (ui)  / F(w)  (7.5.2) 

The  last  term  is  an  error  which  must  be  removed  or  at  least  minimised. 

Here  advantage  may  be  taken  of  the  randomness  of  the  noise  which  implies 
that  its  mean  is  zero.  Hence  if  many  estimates  of  the  "noisy"  transfer 
function  HN(w)  are  taken  and  then  averaged,  the  mean  result  is  a better 
estimate  of  the  true  transfer  function  H(w)  as  the  error  term  should 
average  to  zero.  This  of  course  requires  the  noise  to  be  truly  random 
and  that  a sufficiently  large  number  of  estimates  be  taken. 
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Next  consider  the  effect  of  this  noise  on  the  power  spectrum  of 
the  response.  From  eqn  (7.4.2) 


GXX <w)  = F(u)  + (H(u).  F(u)  + N(u>))* 


|h(o>)  | 2 GFF(a>)  + H (u> ) + H*  (to)  G^fm)  + G^Od) 


(7.5.3) 


I 


where  the  last  three  terms  represent  '.he  error  due  to  the  noise.  The  two 
cross-power  spectra  G_,  (w)  and  G (w)  are  measures  of  the  relationship 
between  the  input  force  F (u>)  and  the  noise  N (u>)  (they  are  in  fact  the 
Fourier  transforms  of  the  cross-correlations  of  these  two  signals) . 

Since  there  is  no  such  relationship  these  two  spectra  have  random  values 
and  will  average  to  zero  if  sufficient  estimates  of  G^lu)  are  taken  and 

averaged.  The  last  term  G (to)  is  a measure  of  the  energy  of  the  noise 

at  each  frequency  and  so  will  not  average  to  zero  but  will  always  remain 
as  an  error  term.  If  we  denote  this  mean  value  of  G^(w)  by  Gxx(w)  we 


5^(0))  = |h(u))|2  gff(o»)  + G^tm) 


(7.5.4) 


Applying  eqn  (7.4.4)  to  obtain  the  modulus  of  the  transfer  function  we 


|HN(u,)|2  = Gxx(o>)  / GpF  ((D) 

= | H (u>)  | 2 + G (u>)  / G (0>) 


(7.5.5) 


The  error  term  here  is  always  positive  so  that  the  estimated  transfer 
function  is  always  larger  than  the  true  one.  However,  if  the  noise  and 
input  force  have  roughly  constant  energy  over  the  frequency  range,  then 
the  function  | («u)  | 2 may  still  have  most  of  the  essential  characteristics 

of  the  true  function  1 H (uj)  | 2 . That  is,  the  peaks  representing  resonances 
may  still  be  apparent,  especially  if  the  darping  is  low  and  qualitative 
estimates  of  this  damping  may  be  made  from  the  "sharpness"  of  the  peaks. 
If,  however,  either  Gm,(u)  or  G (w)  varies  greatly  over  the  frequency 
range,  |hn(u>)  | must  be  treated  carefully. 
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Eqn  (7.4.5)  gives  the  expression  for  the  cross-spectrum  as 


GXF(u)  = (H(u)-  F(u>)  F*(u,) 


= H(w).  G^uj)  + G^w) 


(7.5.6) 


Again  this  quantity  G _(<*>)  should  have  random  values  and  average  to  zero, 

NF 


so  the  mean  cross-spectrum  Gxp (w)  is 


gxf(u,)  = H(uK  gff(u) 


(7.5.7) 


whence 


lyw)  = G^w)  / gff(u) 

= H(oj) 


(7.5.8) 


Hence  this  estimate  of  the  transfer  function  is  unaffected  by  uncorrelated 
noise,  provided  only  that  a sufficiently  large  number  of  averages  has  been 
taken. 


This  then  leads  to  the  coherence  function  which  is  a means  of 


determining  both  the  amount  of  noise  present  and  also  whether  sufficient 
averages  have  been  taken.  The  coherence  function  T (u>)  is  defined  as 


r2(o>)  = |gxf((*>)  1 2 / <Gxx(a))-  gff(u>)) 


XX 


FF 


(7.5.9) 


Making  use  of  eqn  (7.4.5)  we  get 


GXF(u)I  =Gvp(U)-  G£f(W) 


XF XF 

X (w) . F*(u) . X* (M) . F(u>) 


Gxx(u)K  gff(u)) 


(7.5.10) 


Hence  r2(o>)  =*  1 for  all  to  irrespective  of  the  noise  content  of  the  signals. 
If  however,  the  averaged  values  of  the  power  spectra  as  given  by  eqns  (7.5.4) 
and  (7.5.7)  are  used  then  the  coherence  function  becomes 


r2(U,)  = / (Gxx(o».  Gff(o») 

-{|h<»)|?  G^p«u)}  / { ( | H Cto)  |2.  Gff(o>)  + G^  (oj)  > 

={  |h(u>)  I?  Gpp (o>)  } /{  |h(u>)  I2.  Gpptu)  + G^lu)} 

< 1 since  G„. (to)  must  be  positive  for  all  values  of  w. 

NN 

Thus,  as  the  number  of  averages  increase  the  value  of  the  coherence 
2 

function  T (m)  decreases  from  unity  until  a minimum  value  is  reached. 

The  smaller  this  minimum  value,  the  greater  the  level  of  noise 

contaminating  the  system  at  that  frequency.  5 

. 
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CHAPTER  8 ; FOURIER  OPTICS  AND  HOLOGRAPHY 


8.1  Introduction 

Fourier  analysis  is  now  being  widely  applied  in  Optics,  f 

especially  in  connection  with  the  theory  of  image  formation,  optical 
data  processing  and  holography.  In  the  following  pages,  a necessarily 
rather  brief  account  is  given  of  some  of  these  applications.  The  key 
idea  throughout  is  the  use  of  the  Fourier  transform  in  handling 
diffraction  or  wave-front  interference  effects. 

9 

Because  optical  systems  are  generally  analysed  in  two  dimensions, 
the  two-dimensional  form  of  the  Fourier  transform  is  required  here.  This 
can  be  written  as  : 

OP 

G(f  ,f  ) - //  g(x,  y)  exp  f-2iri(f  x + f y)1  dx  dy  (8.1.1) 

x y wfl#  x y j | 

with  the  inverse  relation 

g(x,  y)  = IS  G^Jy)  exp  [2iri(fx  x + f y y)]  dfx  dfy  (8.1.2) 

The  physical  significance  of  the  transform  (or  frequency)  parameters 

f , f will  be  taken  up  shortly, 
x y 

While  all  the  formulae  derived  in  the  following  pages  will  be 
cited  in  the  integral  form,  in  practice  the  DFT  representation  is  generally 
used  for  actual  calculations. 

It  will  be  necessary  throughout  to  assume  several  results  from  Optics 
theory  with  relatively  little  discussion.  Details  of  these  may  be  found 
in  standard  texts  e.g.  those  of  Goodman  (1),  Shulman  (2)  and  Fowles  (3). 

8.2  Frequency  Analysis  of  Optical  Systems  i 


In  this  discussion  only  coherent  optical  systems  will  be  considered 
i.e.  those  systems  illuminated  by  light  which  has  the  property  of  being 
able  to  form  interference  fringes. 


The  term  "frequency"  used  in  the  following  sections  does  not 
refer  to  the  phase  variations  with  time  of  the  light  itself,  but 
refers  to  the  spatial  or  angular  frequency  at  a particular  plane. 
The  various  spatial  frequency  components  in  a light  beam  may  be 
considered  as  plane  waves  propagating  at  different  angles. 


From  the  one-dimensional  view  in  Fig.  8.1(a),  a plane  mono- 
chromatic wave  strikes  a flat  surface  at  an  angle  8 to  the  perpen- 
dicular to  the  surface.  The  lines  of  zero  phase  are  X apart  along 
the  direction  of  propagation,  and  L apart  alon^  the  plane.  The 
spatial  frequency  along  the  plane  is  given  by  : 


Considering  the  two-dimensional  case,  (Fig.  8.1(b)),  if  the  plane 
wave  strikes  the  flat  surface  obliquely,  such  that  the  direction  of 
wave  propagation  makes  angles  with  the  perpendicular  to  the  surface  of 
9 in  the  x direction  and  a in  the  y direction,  then  the  spatial 
frequency  components  cure  given  by 


Consider  a monochromatic,  but  otherwise  unspecified  complex  light 
wave  propagating  in  the  positive  z direction,  as  defined  in  Fig.  8.2 
At  the  x,  y plane,  the  light  field  g(x,  y,  o)  has  a spectrum 


A unit  amplitude  plane  wave,  propagating  WiCfi  direction  cosines  (a,  8,  y) 
is  described  (1)  by 


180  - 


B(x,  y,  z)  = exp  [2vi(ax  + By  + yz)  /xl  (8.3.2) 


Y = A.  - a?  - B2' 

Comparing  the  exponential  terms  of  (8.3.1)  and  (8.3.2)  it  might  seem 
reasonable  that  the  light  field  g(x,  y,  o)  could  be  regarded  as  a plane 
wave  with  direction  cosines 

a = Xf  , B = Xf 

x y 

and 

y = / {1  - (Xfx)2  - (Xfy)2} 

This  is  indeed  borne  out  by  a more  exact  analysis  which  takes  as  its 
starting  point  the  Kirchhoff  diffraction  formula  (4) . 

Eqn  (8.2.1)  may  then  be  re-written  as 
00 

G (a/X , B/X)  = IS  g(x,  y,  o)  expjj-2iri  (ax/X  + By/X  )^j.  dxdy 
—00 

(8.3.3) 

This  is  known  as  the  angular  or  spatial  frequency  spectrum  of 
g(x,  y,  o). 

Now  assume  an  opaque  screen  in  the  x,  y plane  with  an  aperture  E 
(Fig.  8.2)  which  is  large  compared  with  X.  The  transmittance  of  E is 
defined  as 

t(x,  y)  = 1 within  E 
t(x,  y)  =0  elsewhere. 

Applying  the  classical  Kirchhoff  boundary  conditions  (4) 

(i)  Within  E,  the  incident  beam  is  the  same  as  it  would  be  without 
the  screen  being  present. 


The  beam  immediately  behind  the  screen  is  given  by 


gt(x,y,o)  = gi(x,y,o)  t(x,y) 


(8.3.4) 


where  is  the  transmitted  beam  and  g^  is  the  incident  beam. 


Since  a product  in  the  space  (x,y)  domain  becomes  a convolution 
in  the  spatial  frequency  domain,  (8.3.4)  may  be  written  in  terms  of 
spatial  frequency 


Gt(o/A,8/X)  = Gi(a/A,  B/A)  * T (a/A,  B/A)  (8.3.5) 


where 

00 

T (a/A , B/A)  = If  t(x,y)exp{-2iri(ax/A  + By/A)}  dxdy  (8.3.6) 

— oo 


The  diffracted  spectrum  is  thus  found  by  convolving  the 
incident  spectrum  with  the  aperture  "spectrum" . The  diffracted  beam 
is  found  by  transforming  the  result  back  into  the  space  domain. 


A thin  bi-convex  lens  may  be  considered  to  perform  a phase 
transformation  on  an  input  light  wave.  The  form  of  this  transformation 
may  be  readily  calculated  from  the  thickness  variations  of  the  lens. 

Using  this  phase  function,  it  can  be  shown  that  a normally  incident  plane 
wave  is  mapped  by  a bi-convex  lens  into  a converging  spherical  wave, 
provided  the  paraxial  approximation  is  valid.  Using  this  result,  it  cm 
be  further  shown,  e.g.  ref.(l),  that  if  an  object  of  transmittance 
tQ(xo,yo),  placed  at  the  front  focal  plane  of  a bi-convex  lens,  is 
illuminated  by  a normally  incident  monochromatic  plane  wave  (Fig.  8.3), 
then  the  lens  output  converges  to  a distribution  at  the  focal  plane  of 
the  lens  which  is  described  by 


is  zero. 


« 


(8.4.1) 


G(xi,yi) 


If 

—00 


Wy.)expl: 


o 1 


ox  u 


where  f is  the  focal  length  of  the  lens  and  the  Fresnel  (near-field) 
approximation  is  assumed.  Note  that  the  linear  co-ordinates  are  used 
for  the  transform,  as  these  are  generally  more  meaningful  to  the  user. 

A generalized  optical  imaging  system  may  be  characterized  by 
the  relation  (1) 


Gi<*i  Vi>  - !S  h<Vy.;  Vy0)  g<J  (*„,*„>  4*0  aY<> 


(8.4.2) 


where  i and  o are  the  image  and  object  subscripts  and  h is  the 
so-called  "impulse  response"  (also  known  as  the  point  response  or 
imaging  kernel),  i.e.  the  response  of  the  system  at  the  image  plane  to 
a point  source  at  the  object  plane. 


A system  may  be  completely  characterized  by  knowledge  of  the 
function  h.  As  an  example,  to  determine  the  value  of  h for  the 
simple  imaging  system  illustrated  in  Fig.  8.3  let  the  object  be  a point 
source,  i.e.  a three  dimensional  6-function  . The  object  produces  a 
spherical  wave,  which,  after  passing  through  the  lens,  is  described  by 

gL(x,y)  = gQ(x,y)  P (x,y)exp  '^2iri  (x2  + y2)/2Xfj  (8.4.3) 

where  P(x,y)  is  the  lens  pupil  function,  defined  as 

p (x,y)  = 1 within  the  lens  aperture 

P(x,y)  = 0 outside  " " " 

and  with  the  paraxial  approximation  for  the  spherical  wave.  The  quantity 
g (x,y)  becomes  the  function  h at  the  image  plane. 

The  complete  expression  for  h at  the  image  plane  is  cumber- 
some, but  with  the  assumption  that  the  image  is  recorded  on  photographic 
film,  which  time  averages  amplitude  to  record  light  intensify  and  thus 
eliminates  phase  terms,  and  defining  M = d . /dQ,  the  expression 
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simplifies  (1)  to 

oo 

h(xi,yi?XQ,yo)  * d/X^d^  /£  P(x,y)  . 

exp  {-2iri  jjx.^  + Mxq)x  + (y^  + MyQ)yJ  /Xd^  dxdy  (8.4.4) 


where  M is  the  system  magnification. 


This  expression  is  the  Fraunhofer  diffraction  of  the  lens 
aperture,  centred  on  the  image  co-ordinates 


x.  = -Mx  ; y.  = -My 
1 o x o 


To  relate  this  result  back  to  the  effect  on  the  image  quality, 
it  is  convenient  to  make  the  following  change  of  variables 


x = -Mx  ; y = -My  ; x = x/Ad . ; y = y/Ad . ; h = h/M 

O O O O 1 1 


and  define 


gg  (xi,yi)  = V"Xi/M'  _yi/M)/M 


(8.4.5) 


as  the  geometric  (no  diffraction)  image  prediction.  The  imaging  process 
is  then  described  by  the  convolution 


9i(xi'yi}  = h(xi'yi)  * VVYi* 


(8.4.6) 


where 


OO 

'^(xi'Yi)  = ^ P (Xd^x,  Xd^y) exp^-2iri (x^x  + y^jJ  dx  dY 


Because  the  function  h(x^,y^)  does  not  have  zero  width,  as 
would  be  the  case  for  a perfect  impulse  response,  the  image  tends  to  be 
smoothed  out,  resulting  in  a loss  of  fine  image  detail. 
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signal  frequency,  are  being  transmitted  at  a rapid  rate.  At  x2»  the 
received  signal  reflected  from  the  scatterer  is 


S2(t)  = B1  exp{2xifR(t-2d/c)}  (8.6.1) 

where  is  an  amplitude  factor  for  the  scatterer  at  x^, 

d = /{d2  + (x2  - xL)2}  * dx  + (x2  - Xx)  2/2dL 

and  c is  the  velocity  of  light. 

With  (x2  - xx)  - va  t and  c = fRXR,  the  return  from  N 
scatterers  at  position  x^  on  the  ground  to  the  aircraft  at  x2  is 
given  by 

S2(t)  = Z Br  exp{2irifR(t-2d1/XRfR  - vat/d1xRfR>  } (8.6.2) 

N 

This  signal  is  synchronously  demodulated  and  the  result  used  to  modulate 
a CRT  trace  to  give  the  range  record.  The  recording  film  is  drawn  past 
the  CRT  face  to  give  the  azimuth  record  (Fig.  8.10).  The  developed  film 
is  read  out  by  being  illuminated  with  a plane  coherent  wave  with 
corrective  optics  to  compensate  for  the  slanted  view  of  the  antenna  and 
the  different  foci  in  range  and  azimuth.  The  fact  that  the  film  record 
in  the  azimuth  direction  is  a phase  and  amplitude  record  formed  by 
comparing  the  reflected  signal  with  a mutually  coherent  reference  oscillator 
on  the  aircraft  means  that  synthetic  aperture  radar  is  closely  related  to 
holography,  to  be  discussed  in  the  next  section. 

All  of  the  principles  discussed  above  apply  equally  well  to  the 
case  of  synthetic  aperture  sonar  imaging. 

8.7  Holography 

8.7.1  Offset  Reference  Holography 

Holography  is  a method  of  imaging  in  which  the  coherent 
light  scattered  by  an  object  is  recorded  interfercmetrically  on  film. 
Illumination  of  the  developed  film  by  coherent  light  yields  an  image  of 


the  object. 


* 


l_ 


I 

! 
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The  most  common  technique  of  formation  is  the  offset 
reference  method  (9),  in  which  a reference  beam  (Fig.  8.11)  is 
directed  onto  a photographic  plate  at  an  angle,  while  the  light  scattered 
(by  reflection  from,  or  transmission  through)  the  object  also  falls  on 
the  plate.  A photographic  record  of  the  interference  fringes  formed  by 
the  two  beams  is  the  hologram. 

Reconstruction  is  achieved  by  illuminating  the  developed 
hologram  with  the  constructing  beam.  Two  images,  one  real  and  one  virtual 
are  obtained  (Fig.  8.12). 

The  parallel  between  the  offset  reference  hologram  and 
Vander  Lugt  filter  formation  (Fig.  8.7)  is  evident.  The  two  concepts 
originated  from  the  same  laboratory  at  about  the  same  time. 

In  the  spatial  frequency  domain,  if  the  object  spectrum 
(Fig. 8. 13 (a) ) is  of  width  w,  then  the  condition  for  separation  of  the 
images  from  the  central  order  beam  (Fig.  8.13(b))  is  that 

a > 3w  (8.7.1) 

where  the  spatial  frequency  of  the  reference,  a = sin  0/X,  where  9 is 
the  reference  beam  angle,  and  X is  the  wavelength. 

The  offset  reference  may  be  considered  an  analogue  of 
the  carrier  in  a radio  transmitter,  with  the  twin  images  being  equivalent 
to  the  two  side  bands. 

8.7.2  Fourier  Transform  Holography 

In  this  method  of  formation,  (10),  am  object,  which 

may  be  considered  as  a collection  of  point  scatterers,  is  illuminated 

by  a coherent  source  (Fig.  8.14).  A mutually  coherent  point  source,  R, 

is  located  a short  distance,  £ , away  from  the  centre  of  the  object,  and 

o 

coplanar  with  it.  The  spherical  wave  at  the  recording  ^ilm  due  to  the 
reference  source  is  approximated  by 

Ar(x)  * Aj^  exp (2iri  x2/2Xf)  (8.7.2) 
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The  wave  at  the  film  plane  from  an  object  point  at  is  given  by 


Aq(x)  = A(So>  exp{2iri(x  - £q)2/  2XfJ 


(8.7.3) 


If  5 is  small  so  K « x,  then  the  whole  object  wave  is  described 
o o 


AQbj  (x)  = J A(£)  expjjhri  x2/2XfJexp[-2iri  x£/fxj  d£  (8.7.4) 


The  intensity  recorded  by  the  film  is 


K*>  - (*R  + W (ar  + AObj>* 


| All  2 + |Aobj|2  + A XJ  A(C)  exp(-2iri  x£/fX)d£ 

?1 


+ A.  / A(£)*  exp(+2xi  xS/Xf)d£ 


(8.7.5) 


where  the  asterisk  denotes  the  conjugate.  The  last  two  terms  of  (8.7.5) 
are  the  Fourier  transform  of  the  object  and  its  conjugate;  hence  the 
name  Fourier  Transform  holography. 


Reconstruction  is  achieved  by  illuminating  the  hologram 
with  the  coherent  source  and  converging  the  output  with  a lens  (Fig.  8.15) 
to  produce  twin  offset  images. 

The  Fourier  transform  holography  is  restricted  to  two- 
dimensional  objects,  but  has  a relatively  low  spatial  frequency  content, 
especially  if  is  small,  enabling  normal  film  to  be  used.  The  offset 

reference  type  of  hologram,  in  contrast,  can  produce  three-dimensional 
images  but  the  high  spatial  frequency  content  requires  spectroscopic  plate 
for  recording. 


' 
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FIG.  8.3  GEOMETRY  FOR  AN  EXACT  FOURIER  TRANSFORM  OF  AN  OBJECT 
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FIG.  8.4  GEOMETRY  TO  FOURIER  TRANSFORM  MESH  AND  INVERSE 
FOURIER  TRANSFORM  ITS  SPECTRUM 
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PIG.  8.5  FOURIER  TRANSFORM  OF  MESH 


FIG.  8.6  COHERENT  PROCESSOR  GEOMETRY 
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FIG.  8.8  OUTPUT  PLANE  INTENSITY  OF  THE  COHERENT  PROCESSOR 
WITH  VANDER  LUGT  FILTER 
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FIG.  8.12  OFFSET  REFERENCE  HOLOGRAM  RECONSTRUCTION 
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